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PREFACE 


The  work  reported  herein  was  sponsored  by  the  Office,  Chief  of 
Engineers  (OCE) ,  US  Army  through  the  US  Army  Corps  of  Engineers  Waterways 
Experiment  Station  (WES)  as  a  pilot  study  in  the  reevaluation  and  possible 
updating  of  the  system  currently  used  for  predicting  soil  traff icability 
through  the  estimation  of  soil  moisture  content.  The  study  was  performed 
July  1981  through  July  1982  and  was  authorized  by  WES  Work  Order  No.  DACty- 
39-81>*M— 3793,  dated  13  July  1981.  Funding  for  the  study  was  provided  by 
Research  and  Development  Project  No.  4A161102AT24,  Task  Area  E3,  Work 
Unit  004,  "Weather  and  Climate  Influence  on  Mobility."  The  work  was 
performed  by  the  Department  of  Soil  Science  and  Biometerology,  Agricultural 
Experiment  Station,  Utah  State  University;  Dr.  L.  F.  Hall  was  the  Principal 
Investigator. 

The  work  was  monitored  by  Mr.  C.  J.  Nuttall,  Jr.,  Chief,  Mobility 
Systems  Division  (MSD) ,  Geotechnical  Laboratory  (GL) ,  WES,  and  by  Mr.  G.  W. 
Turnage,  MSD,  under  the  general  supervision  of  Dr.  W.  F.  Marcuson  III, 
Chief,  GL. 

COL  Tilford  C.  Creel,  CE,  and  COL  Robert  C.  Lee,  CE,  were  Commanders 
and  Directors  of  WES  during  the  conduct  of  this  study  and  preparation  of 
this  report.  COL  Creel  was  also  the  Contracting  Officer.  Mr.  Fred  R. 

Brown  and  Dr.  Robert  W.  Whalin  were  Technical  Directors. 
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A  REVIEW  OF  DEVELOPMENTS  IN  AGRICULTURAL  SCIENCE  APPLICABLE 
TO  MILITARY  SOIL  MOISTURE  PREDICTION  REQUIREMENTS 


PART  Is  INTRODUCTION 


Background 


1.  Soil  moisture  content  is  critical  to  military  hydrology  and 
mobility.  Excessive  soil  moisture  reduces  infiltration  of  rain  or  snowmelt 
water  and  decreases  soil  strength.  These  effects  lead  to  high  stream  levels, 
flooding,  reduced  soil  traff icability,  and  general  impairment  of  activity. 
Prediction  of  the  moisture  content  of  soils  is  therefore  important  to 
effective  planning  of  field  operations. 

2.  The  soil  moisture  prediction  methods  that  are  currently  used  by 
the  US  Army  for  hydrology  and  mobility  models  are  based  on  work  performed  in 
the  1930's  and  1950's,  respectively.  While  age  is  not  of  itself  a  condemna¬ 
tion  of  a  model,  witness  the  Pythagorean  theorem,  and  while  some  improvements 
in  the  models  have  been  made  since  their  formulation,  the  present  soil 
moisture  models  for  military  applications  are  distinctly  inferior  to  the 
state  of  the  art. 

3.  Soil  moisture  prediction  models  for  entire  watersheds  and 
individual  fields  have  been  developed  in  the  civilian  community  since  the 
current  military  models  were  formulated.  The  more  advanced  models  directly 
solve  the  differential  equations  related  to  soil  water  movement  by  numerical 
integration  with  digital  computers.  These  newer  physically  based  models, 
therefore,  are  in  dramatic  contrast  with  the  correlation  methods  used  in  the 
military  models.  The  newer  models  are  possible  because  of  overwhelming 
increases  in  the  computing  power  of  modern  machines,  compared  to  past 
decades,  and  because  of  fundamental  advances  in  evapotranspiration  and  soil 
physics  theory.  Modern  models  have  also  benefited  from  improvements  in 
computational  methods,  and  from  both  improved  field  measurement  techniques 
and  field  data  sets. 

4.  Recent  reviews  of  the  military  soil  moisture  prediction  models 
have  noted  the  wealth  of  knowledge  which  exists  in  the  civilian  community  on 
this  subject,  e.g..  Proceedings  of  the  Military  Hydrology  Workshop,  17-19  May 


1978,  Vicksburg,  Mississippi  (US  Army  Engineer  Waterways  Experiment  Station 
(WES),  1979).  Evaluation  of  civilian  methods,  improved  technology  transfer, 
and  consideration  of  military  model  updating  have  generally  been  recommended 
by  the  review  panels.  Because  much  of  the  recent  modeling  work  has  been 
performed  in  connection  with  agricultural  problems,  such  as  irrigation 
scheduling,  yield  prediction,  and  runoff  estimation,  the  present  review  is  a 
large  step  toward  those  goals. 


Objectives 

5.  The  objectives  of  the  research  discussed  herein  were  to  review 
the  soil  moisture-weather-time  relations  developed  in  agricultural  science, 
and  to  report  on  the  applicability  of  such  relations  to  military  traffic- 
ability  problems.  Emphasis  was  to  be  placed  on  the  effects  of  soil  water — 
its  movement,  replenishment,  and  depletion  by  evapotranspiration. 


Scope 


6.  The  project  began  with  a  thorough  review  of  the  current  soil 
moisture  prediction  model  used  for  military  traff icability  problems,  and  the 
historical  development  of  the  model.  The  review  included  interviews  with 
WES  personnel  familiar  with  the  model  and  its  history,  and  a  review  of  WES 

documents  written  by  both  WES  and  US  Forest  Service  personnel  during  model 

\ 

development  (see  Meyer,  1976,  for  Bibliography). 

7.  It  became  clear  early  in  the  project  that  numerical  modeling  of 
soil  moisture  was  the  activity  which  would  provide  the  best  route  through 
the  literature,  and  was  the  most  significant  development  in  the  last  three 
decades.  Other  theoretical  and  experimental  developments  were  identified 
and  reviewed  as  they  related  to  modeling.  Several  approaches  to  the  litera¬ 
ture  were  included  in  a  concerted  effort  to  ensure  that  each  major  modeling 
group  was  discovered  and  represented: 

a.  A  review  paper  on  crop  yield  models  (Kanemasu  et  al.,  1980) 
and  a  review  paper  on  soil  moisture  determination  methods 
(Schmugge,  Jackson,  and  McKim,  1979)  were  scanned  for 
appropriate  references.  The  latter  included  a  review  paper 
on  soil  moisture  modeling  (Hildreth,  1978)  that  provided 
additional  references.  (A  journal  paper  version  of  the  second 
reference  (Schmugge,  Jackson,  and  McKim,  1980)  was  discovered 
during  subsequent  searching.) 

A 


_b.  Recent  indices  for  the  Soil  Science  Society  of  America  Journal 
(formerly  Proceedings) ,  Soil  Science,  Transactions  of  the 
American  Society  of  Agricultural  Engineers,  and  Soviet  Soil 
Science  were  reviewed  under  title  and  subject  categories  for 
relevant  publications. 

£.  The  table  of  contents  for  each  journal  from  which  a  paper  was 
obtained  was  read  in  search  of  relevant  titles  or  work  by 
authors  known  to  be  involved  in  soil  moisture  modeling. 

Several  papers  not  otherwise  referenced  were  thereby  discovered. 
While  not  comprehensive,  this  technique  included  most  of  the 
recent  work  (to  December  1980)  published  in  Soil  Science, 

Soil  Science  Society  of  America  Journal,  Water  Resources 
Research,  and  Transactions  of  the  American  Society  of 
Agricultural  Engineers. 

(1.  Relevant  referenced  books  which  were  available  in  the  WES 
library  were  reviewed. 

£.  The  Defense  Technical  Information  Center  literature  search 
capability  was  accessed  via  a  library  computer  terminal 
through  a  number  of  key  words  related  to  soil  moisture 
modeling. 

8.  The  scope  of  this  literature  review  was  therefore  large,  but  not 
exhaustive.  Work  published  since  December  1980  was  generally  not  included. 
Modeling  effort  which  has  not  been  revealed  through  publications  or  reference 
in  the  materials  searched  may  well  exist.  Also,  some  unobtainable  work  has 
been  noted,  particularly  work  of  the  modeling  group  at  Wageningen,  The 
Netherlands.  These  caveats  notwithstanding,  the  literature  search  conducted 
was  fully  adequate  to  the  purpose  of  the  project  in  that  it  has  identified  a 
number  of  approaches  to  soil  moisture  modeling  which  establish  the  applica¬ 
bility  of  these  developments  to  military  problems. 

9.  The  materials  identified  as  relevant  to  this  project  included  the 


following: 


All  references  to  numerical  modeling  of  soil  moisture,  water 
infiltration,  redistribution,  drainage,  and  evapotranspira- 
tion. 

All  references  to  numerical  modeling  of  crop  yields  as  due  to 
soil  moisture  availability  that  also  included  a  significant 
treatment  of  soil  moisture. 

References  to  analytic  solutions  of  soil  moisture  problems 
that  were  referenced  as  numerical  modeling  check  values. 


d_.  References  to  procedures  for  the  calculation  of  hydraulic 
conductivity  or  moisture  diffusivity  from  particle  size 
distributions,  pore  size  distributions,  or  the  soil  charac¬ 
teristics  curve  of  moisture  content  versus  matric  potential 
(soil  moisture  tension) . 


e.  References  to  characteristic  curves  that  are  typical  for 
soils  identified  by  type  (in  contrast  to  specific  samples). 

f.  Principal  references  to  measurement  and  estimation  of  soil 
property  heterogeneity  within  single  fields  or  soil  types. 

£.  Survey  papers  on  soil  moisture  determination  via  remote 
sensing. 

h.  References  identifying  field  data  that  may  be  used  for  model 
validation  or  that  have  been  so  used. 

i.  References  to  soil  moisture  modeling  with  respect  to  soil 
survey  maps,  and  large  area  modeling. 

Principal  references  to  modeling  of  water  tables,  ground 
water,  and  saturated  flow. 

10.  These  reference  materials  have  been  reviewed  under  the  headings, 
"Soil  Physics,"  "Evapotranspiration, "  and  "Numerical  Modeling  of  Soil  Moisture" 
in  Appendices  A,  B,  and  C,  as  well  as  being  utilized  in  preparation  of  the 
body  of  this  report. 

Definitions 


11.  Several  terms  used  in  this  report  are  given  here  with  general 

definitions  for  the  convenie  ce  of  the  reader.  Appropriate  texts  should  be 

consulted  for  more  complete  and  precise  definitions  of  technical  terms. 

Accretion.  Addition  of  moisture  to  a  soil  column. 

Algorithm.  A  set  of  steps  to  be  performed  to  solve  a  problem, 
used  for  mathematical  steps  in  a  numerical  solution. 

Analytic  solution.  An  exact  mathematical  expression  for  the 
state  of  a  system  throughout  space  and  time.  It  is  bounded  in 
accuracy  by  the  validity  of  assumptions  needed  to  derive  it  and 
the  accuracy  of  input  variables  and  parameters  used  in  the 
expression. 

Catenary  process.  A  total  process  made  up  of  a  series  of  partial 
processes  in  which  the  total  process  rate  is  limited  by  the 
slowest  partial  process  rate. 

Correlation  relations.  Empirical  solutions  derived  by  regression 
of  assumed  dependent  variable  data  against  assumed  independent 
variable  data. 

Crop  coefficient.  A  multiplier  of  calculated  potential  evapotrans¬ 
piration  used  to  account  for  crop  characteristics  which  differ 
from  short  grass,  including  crop  height  and  stage  of  growth. 

Depletion.  Removal  of  moisture  from  a  soil  column. 


Desorption.  The  movement  of  water  out  of  a  soil  sample  due  to 
soil  moisture  tension  (matric  potential) . 

Drainage.  The  movement  of  water  out  of  a  soil  column  due  to 
gravity. 

Drying  curve.  The  soil  characteristic  curve  derived  during 
drying  of  a  saturated  soil  sample. 

Empirical  solution.  An  approximate  expression  derived  from 
experimental  data  for  the  relationship  between  variables. 

Evaporation.  The  loss  of  soil  moisture  or  ponded  water  at  or 
through  the  soil  or  water  surface  after  vaporization,  including 
vaporization  of  dew  or  intercepted  precipitation  from  plant 
surfaces. 

Evapotranspiration.  Evaporation  plus  transpiration. 

Explicit  numerical  method.  Dependent  variables  at  a  given  time 
step  of  a  numerical  solution  are  calculated  from  known  values  of 
dependent  variables  that  were  calculated  at  a  prior  time  step. 

Field  capacity.  The  moisture  coutent  which  prevails  in  a  soil 
column  after  a  period  of  drainage  from  saturated  conditions, 
nominally  two  days  with  no  evapotranspiration. 

Field  heterogeneity.  Variability  in  space  of  a  soil  property. 

Hydraulic  conductivity.  The  ratio  of  water  flux  density  to 
hydraulic  potential  gradient. 

Hydraulic  head.  See  hydraulic  potential. 

Hydraulic  potential.  Sum  of  matric  and  gravitational  potentials 
in  unsaturated  systems. 

Hysteresis.  A  system  property  for  which  the  state  of  a  system 
and  changes  to  that  state  depend  on  the  direction  of  change  and 
on  previous  states  of  the  system. 

Implicit  numerical  method.  Dependent  variables  at  a  given  time 
step  of  a  numerical  solution  are  calculated  simultaneously  with 
other  dependent  variable  values  at  other  spatial  points  for  the 
same  time  step.  The  algorithm  may  include  known  values  of 
dependent  variables  from  prior  time  steps,  also. 

Infiltration.  The  entry  of  water  into  a  soil  column  through  the 
soil  surface. 

Integro-dif ferential  equation.  An  equation  in  which  terms  of  both 
differential  and  integral  form  appear. 

Latent  evaporation.  The  evaporation  rate  from  an  evaporimeter , 
an  instrument  with  a  saturated  porous  surface  which  is  exposed  to 
sun  and  wind. 

Lidar .  Light  detection  and  _ranging.  Active  remote  sensing  system 
using  laser  illumination  and  telescopic  receiver  to  measure 
distance  to  and  reflection  by  dust,  clouds,  or  precipitation 
events. 
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Matric  potential.  The  pressure  (negative)  required  to  hold  water 
in  equilibrium  with  soil  moisture  across  a  porous  membrane  due  to 
adsorptive  forces  on  the  soil  water.  (See  soil  physics  text  for 
details . ) 

Moisture  diffusivity.  The  ratio  of  water  flux  density  to  the 
gradient  of  soil  moisture  content. 

Nonlinear .  An  adjective  used  to  denote  systems  and  their 
mathematical  descriptions  which  involve  products  of  variables  and 
their  derivatives  or  other  functions  of  those  variables,  e.g., 
when  conductivity  is  a  function  of  potential,  then  the  product  of 
conductivity  and  the  gradient  of  potential  is  nonlinear. 

Numerical  solution.  An  approximate  algorithm  through  which  the 
state  of  a  system  may  be  calculated  for  discrete  points  in  space 
and  time  via  stepwise  representation  of  temporal  and  spatial 
gradients.  Normally,  compared  to  analytic  solutions,  numerical 
solutions  are  less  bounded  by  system  simplifications,  are  equally 
bounded  by  input  variable  and  parameter  accuracy,  and  are  more 
bounded  by  computational  expense  required  for  stability  and 
equivalent  accuracy. 

Permanent  wilting  point.  The  moisture  content  which  prevails  in 
a  column  of  soil  when  plants  are  observed  to  wilt  permanently. 

Porosity.  The  ratio  of  total  pore  volume  to  bulk  soil  volume. 

Potential  evapotranspiration.  The  rate  of  evapotranspiration 
that  would  occur  at  a  well  watered  surface,  i.e.,  limited  only  by 
solar,  atmospheric,  and  submedium  influences  which  control  the 
supply  of  latent  heat  for  vaporization  and  remove  water  vapor. 

Pressure  head.  See  matric  potential. 

Redistribution.  The  movement  of  soil  water  within  a  soil  column 
in  response  to  hydraulic  gradients,  normally  applied  to  the 
period  following  infiltration. 

Remote  sensing.  Noncontact  measurement  using  properties  of 
electromagnetic  radiation  and  its  emission  or  reflection  by 
surfaces. 

Scanning  curve.  One  of  an  unbounded  set  of  soil  characteristic 
curves  derived  during  rewetting  of  a  partially  dried  soil  sample, 
or  drying  of  an  initially  unsaturated  soil  sample.  Scanning 
curves  run  between  wetting  and  drying  curves  of  the  sample. 

Soil  characteristic  curve.  Graph  of  matric  potential  against 
soil  moisture  content. 

Soil  moisture  content.  Water  present  in  the  soil  matrix.  May  be 
expressed  as  volume  of  water  per  volume  of  soil,  equivalent 
depths  of  water  and  soil,  or  mass  of  water  per  mass  of  soil. 

Soil  moisture  model.  A  set  of  relationships,  normally  mathematical 
expressions,  describing  soil  moisture  content  in  space,  and  time 
as  a  consequence  of  water  input,  water  movement,  water  storage, 
and  water  loss.  The  solution  technique  is  considered  part  of  the 
model,  particularly  in  the  case  of  numerical  solutions. 
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considered  constant  during  calculation  of  matric  potential  for  one  time 
step,  a  new  hydraulic  conductivity  is  computed  with  the  calculated  matric 
potential  for  the  time  step,  and  the  loop  is  continued  until  two  values  of 
matirc  potential  for  the  time  step  are  sufficiently  close.  Calculation  then 
proceeds  to  the  next  time  step  after  all  grid  points  have  been  resolved. 

This  method  has  been  used,  though  not  extensively,  in  soil  moisture  models. 

49.  The  more  common  approach  to  the  problem  of  nonlinearity  is  to 
remove  the  nonlinearity  by  use  of  an  estimated  hydraulic  conductivity  that 
is  considered  constant  for  the  time  step.  This  may  be  done  explicitly  or 
implicitly.  Explicit  linearization  of  the  equation  uses  conductivity  (or 
diffusivity  or  specific  water  capacity,  depending  on  the  formulation  of  the 
basic  equation — see  Appendix  A)  that  is  calculated  via  known  values  of 
matric  potential  or  soil  moisture  content.  Usually,  a  spatial  average  value 
across  one  or  two  grid  intervals  is  used  with  improved  results  in  comparison 
to  a  single  point  value.  Alternately,  a  predicted  value  of  matric  potential 
for  one-half  step  forward  in  time  can  be  calculated  using  an  explicitly 
linearized  equation  (predictor  step).  This  predicted  value  is  used  to 
obtain  the  half-step  value  of  hydraulic  conductivity,  which  is  then  considered 
constant  for  the  whole  step  and  used  to  calculate  a  corrected  value  of 
matric  potential  for  a  whole  time  step  (corrector  step) .  In  addition  to 

this  predictor-corrector  approach,  a  number  of  other  estimates  of  the  matric 
potential  after  a  time  step  have  been  used  by  various  investigators  to 
calculate  the  coefficients  for  use  over  the  time  step.  The  method  used  to 
deal  with  nonlinearity  should  be  noted  for  each  model. 

50.  Numerical  models  may  use  a  method  of  removing  the  nonlinearity  of 
the  equations  because  each  time  step  of  the  calculations  is  essentially 
independent.  The  results  oi  prior  steps  may  be  considered  initial  conditions 
for  a  discrete  change.  Analytic  models,  on  the  other  hand,  treat  continuous 
changes  via  functional  relationships  which  do  not  permit  introduction  of 
temporarily  constant  coefficients.  The  discrete  stepping  of  numerical 
solutions  is  therefore  of  distinct  benefit  in  treating  nonlinear  systems,  in 
compensation  for  less  accurate  representation  of  continuous  changes.  Also, 
hysteresis  is  easily  treated  by  numerical  models,  because  history  cau  be 
considered  for  a  step  and  the  proper  values  of  conductivity,  matric  potential, 
etc.,  can  be  selected  from  functional  relations  or  tables.  Analytic  models 
(and,  generally,  empirical  models)  do  not  possess  this  flexibility. 
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as  central  versus  one-sided  differencing,  handling  of  initial  and  boundary 
conditions,  and  discussion  of  theoretical  accuracy  of  various  methods  of 
differencing,  are  not  treated  here.  The  unfamiliar  reader  is  referred  to 
Remson,  Hornberger,  and  Molz  (1971),  or  to  general  texts  on  numerical  methods. 

45.  Explicit  finite  differencing  results  in  an  equation  for  the  value 
of  a  dependent  variable  at  a  succeeding  time  step  and  spatial  grid  point  in 
terms  of  known  values  of  dependent  and  independent  variables  and  coefficients. 
The  equation  can  be  directly  solved  for  the  dependent  variable.  Explicit 
finite  differencing  is  the  most  straightforward  approach  and  the  simplest  to 
program  and  may  be  used  to  resolve  the  entire  spatial  grid  one  point  at  a 
time.  However,  small  time  steps  must  be  used  based  on  spatial  grid  size 

and  process  physics,  to  avoid  instability. 

46.  Implicit  finite  differencing  results  in  a  system  of  equations  for 
the  values  of  a  dependent  variable  at  all  spatial  grid  points  at  a  succeeding 
time  step  in  terms  of  the  dependent  variable  at  adjacent  grid  points.  The 
system  of  equations  results  because  spatial  gradients  of  the  dependent 
variable  are  written  using  future  values  of  the  variable.  While  not  as 
simple  to  program  or  visualize  as  an  explicit  method,  an  implicit  method  has 
advantages.  It  is  more  stable  and  more  economical  for  a  single  time  step 
resolution  of  the  entire  spatial  grid,  and  its  stability  permits  further 
economy  through  larger  permissible  time  steps. 

47.  One  implicit  finite  difference  approach  that  has  been  used 
extensively  is  the  Crank-Nicolson  scheme.  In  this  scheme,  spatial  gradients 
are  written  as  averages  of  gradient  terms  involving  known  and  unknown  values 
of  the  dependent  variable,  i.e.,  explicit  and  implicit  terms.  While  more 
complex  in  formulation  than  a  fully  implicit  form,  the  Crank-Nicolson  scheme 
results  in  improved  accuracy  due  to  better  control  of  problems  which  arise 
within  the  numerical  integration  itself.  The  integration  over  a  time  step 
is  performed  using  a  gradient  appropriate  to  the  step  midpoint  due  to  the 
averaging  of  explicit  and  implicit  terms. 

48.  The  Richards  equation  (and  each  of  its  equivalent  forms)  is  a 
strongly  nonlinear  partial  differential  equation.  The  nonlinearity  arises 
since  Darcy's  law  is  applied  to  uusaturated  media  for  soil  moisture  predic¬ 
tion,  and  hydraulic  conductivity  (coefficient)  is  a  strong  function  of 
matric  potential  (gradient  term)  in  unsaturated  soils.  The  nonlinear  equation 
may  be  solved  by  iteration,  where  an  estimated  hydraulic  conductivity  is 


42.  Matric  potential  and  volumetric  moisture  content  are  both  unknowns 
for  models  of  the  unsaturated  zone,  and  the  Richards  equation  contains  both 
(see  Appendix  C,  paragraph  6).  One  of  the  two  is  normally  eliminated  from 
the  equation  mathematically,  and  then  obtained  from  the  soil  characteristic 
curve  after  the  other  has  been  calculated  by  the  numerical  solution.  Elimina¬ 
tion  of  volumetric  moisture  content  enables  more  general  application  of  the 
model,  as  it  is  then  applicable  to  cases  with  ponded  water  on  the  surface. 
Further,  when  total  potential  is  used  as  the  unknown,  the  equation  also 
applies  to  saturated  flow.  On  the  other  hand,  elimination  of  matric  potential 
leads  to  a  diffusion  equation  with  advantages  in  some  applications.  Both 
approaches  have  been  used  with  good  results  for  appropriate  initial  and 
boundary  conditions,  although  for  general  modeling,  it  is  better  to  use  the 
matric  potential  form  of  the  equation  (Philip,  1958). 

43.  Given  the  Richards  equation  (or  the  equivalent  in  either  matric 
potential  or  volumetric  moisture  content)  as  the  basic  mathematical  descrip¬ 
tion  of  the  physical  process  of  soil  moisture  movement  and  storage,  or  the 
principle  of  conservation  with  Darcy's  law  flow  into  and  out  of  soil  volumes, 
the  equations  must  be  rewritten  in  a  form  suitable  for  computer  solution. 

The  most  commonly  used  approach  is  to  rewrite  gradients  in  space  and  time  in 
terms  of  finite  spatial  and  temporal  steps,  i.e.,  as  finite  differences. 

The  solution  is  then  calculated  for  the  points  of  the  discrete  spatial  grid 
by  stepping  through  time.  Step  sizes  must  be  selected  to  achieve  an  optimum 
balance  between  solution  accuracy  and  computational  expense  and  time  for 
each  situation,  with  additional  consideration  required  to  ensure  that  the 
solution  is  stable  and  converges  to  the  "true"  result.  An  alternate  method 
coming  into  more  frequent  use  is  to  use  principles  of  variational  calculus 
to  formulate  the  numerical  algorithm,  frequently  resolving  the  spatial 
variation  in  terms  of  finite  elements  rather  than  at  points  of  a  grid.  This 
approach  is  well  covered  in  Remson,  Hornberger,  and  Molz  (1971)  and  Guymon, 
Scott,  and  Herrmann  (1970).  At  this  time,  it  appears  most  suitable  to 
areas  with  complex  boundaries.  Because  most  agricultural  applications  have 
involved  only  the  finite  difference  approach,  it  has  been  emphasized  in  this 
report. 

44.  Two  other  major  decisions  complete  the  outline  of  a  specific 
model,  namely,  explicit  versus  implicit  finite  differencing  and  the  method 
of  treating  nonlinearity  of  the  equations.  Other  details  of  solutions,  such 
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39.  Judgements  of  various  models  given  in  the  above-cited  references 
must  be  viewed  in  the  context  of  the  respective  texts  or  reports,  however. 

Both  crop  yield  models  and  remote  sensing  techniques  are  too  coarse  to  be 
able  to  use  the  information  from  the  more  sophisticated  models;  thus  they 
tend  to  emphasize  the  economy  of  budget  models  for  soil  moisture  prediction. 
While  runoff  estimation  may  fall  in  a  similar  category  with  respect  to 
resolution  in  depth  and  area  as  yield  models  and  remote  sensing  techniques, 
the  depth  and  areal  resolution  for  soil  strength  prediction  must  be  greater. 
This  is  particulary  true  when  mobility  is  considered  beyond  a  simple  GO-NOGO 
prediction,  or  when  hours  count  in  tactical  operation,  or  when  the  traffic- 
ability  predictions  apply  for  vehicles  covering  a  broad  range  of  traffica- 
bility  capabilities. 

40.  Numerical  models  for  soil  moisture  prediction  are  based  on  Darcy's 
law  (moisture  flow  is  proportional  to  hydraulic  potential  gradient)  and  the 
principle  of  the  conservation  of  mass  (equation  of  continuity) .  Models  have 
been  developed  which  are  based  on  the  numerical  solution  of  a  single  equation 
of  continuity  (see  Appendix  A) ,  and  also  based  on  the  change  of  water  sub¬ 
stance  in  a  soil  layer  (or  volume)  from  the  sum  of  flows  across  layer  surfaces. 
Generally,  the  former  approach  has  been  used  by  modelers  who  did  their  own 
programming  in  FORTRAN,  particularly  modelers  using  an  implicit  numerical 
method  of  solution,  while  the  latter  approach  has  been  used  by  most  modelers 
working  with  CSMP  and  using  explicit  numerical  methods.  Neither  representa¬ 
tion  of  process  physics  is  intrinsically  superior  to  the  other  for  model 
applications. 

41.  Flow  in  the  unsaturated  zone  of  the  soil  is  generally  transient 

in  situations  of  interest,  including  temporary  saturation  during  infiltration. 
Soil  moisture  prediction  models  calculate  flow  and  soil  moisture  content  for 
a  future  time  by  stepwise  temporal  change  from  specified  initial  conditions, 
with  specified  boundary  conditions  on  the  volume  for  which  prediction  is 
being  made.  Flow  in  the  saturated  zone  is  generally  treated  as  steady-state 
movement  in  response  to  specified  boundary  conditions,  although  water  table 
movement  (transient)  is  treated  in  a  few  models.  Some  form  of  iterative 
solution  is  used  in  the  saturated  flow  case,  such  that  successive  approxima¬ 
tions  to  the  flow  converge  to  steady-state.  Models  which  treat  both 
unsaturated  and  saturated  zones  require  coupling  between  basically  transient 
and  steady-state  solutions  for  the  respective  zones. 


36.  In  all,  the  concepts  of  potential  evapotranspiration,  catenary 
processes  for  water  movement  through  plants,  and  stages  of  drying  of  a  soil 
have  aided  understanding  of  the  complex  process  of  evapotranspiration.  Most 
soil  moisture  prediction  models  include  these  concepts  and  their  mathematical 
implementation  to  some  degree  at  present. 

Numerical  Models 

37.  Fleming  (1975)  includes  a  short  summary  of  the  development  of 
computing  machines  in  the  beginning  of  his  book  on  watershed  modeling, 
including  the  dramatic  increases  in  their  computational  speed  and  memory 
capacity  during  the  three  decades  since  the  introduction  of  semiconductor 
devices.  The  increases  in  both  speed  and  memory  have  made  possible  models 
of  physical  phenomena  that  were  unthinkable  in  the  early  1950's,  at  least 
outside  science  fiction.  Parallel  software  development  has  produced  computer 
languages,  such  as  FORTRAN,  and  higher  order  software,  such  as  CSMP  (Continu¬ 
ous  System  Modeling  Program) ,  which  greatly  simplify  communication  with  the 
computers.  Numerous  models  have  been  developed  for  soil  moisture  prediction 
which  utilize  these  increased  computational  capabilities,  ranging  from 
simple  mass  balance  models  similar  to  the  WES  soil  moisture  prediction  model 
to  complex  three-dimensional  watershed  models  with  extensive  consideration 

of  process  physics.  The  general  features  of  these  models  are  discussed 
here,  while  Appendix  C  includes  details  of  the  approaches  used  by  several 
active  modeling  groups. 

38.  Several  references  are  particularly  recommended  for  additional 
study;  for  an  overview  of  soil  moisture  prediction  models  with  emphasis  on 
budget  (mass  balance)  approaches — Schmugge,  Jackson,  and  McKim  (1980);  for 

an  excellent  discussion  of  numerical  techniques  in  hydrological  applications, 
including  soil  moisture  prediction  models — Remson,  Hornberger,  and  Molz  (1971); 
for  general  information  on  CSMP  and  a  discussion  of  its  application  to  a  broad 
range  of  specific  soil  moisture  prediction  problems,  including  a  two-dimen¬ 
sional  watershed  model — Hillel  (1977);  for  a  rigorous  comparison  of  six  numeri¬ 
cal  modeling  approaches  to  both  analytic  solutions  and  laboratory  data  for  the 
se  of  infiltration  into  sand — Havercamp  et  al.  (1977);  for  further  details 
on  a  number  of  budget  models — Hildreth  (1978);  and  for  a  review  of  soil  mois¬ 
ture  modeling  as  it  relates  to  crop  yield  models — Kanemasu  et  al.  (1980). 


the  rate  of  evapotranspiration  itself,  as  it  exceeds  the  rate  at  which  water 
can  be  supplied.  Several  developments  have  recently  aided  organization  of 
this  complex  of  factors  with  resulting  progress  in  evapotranspiration 
prediction. 

33.  The  most  important  concept  with  respect  to  evapotranspiration  was 
put  forth  in  the  mid-1940,s  by  Thornthwaite  et  al.  (1944),  namely,  potential 
evapotranspiration.  Potential  evapotranspiration  is  that  evapotranspiration 
rate  that  is  not  limited  by  moisture  supply,  but  only  by  energy  supply  and 
ventilation.  Thornthwaite  (1948)  was  able  to  formulate  useful  regression 
relations  based  on  mean  temperature  and  day  length  for  the  calculation  of 
potential  evapotranspiration,  while  Penman  (1956)  and  Van  Bavel  (1966) 
modeled  potential  evapotranspiration  in  terms  of  the  energy  balance  at  the 
surface  and  ventilation.  This  concept  has  proven  of  value  in  the  separation 
of  limits  to  evapotranspiration  and  more  detailed  study  of  the  astmospheric 
contributions.  Soil  and  plant  limits  to  evapotranspiration  are  then  seen  in 
terms  of  the  difference  between  potential  and  actual  rates. 

34.  A  second  organizing  concept  of  great  importance  was  put  forth  by 
Van  den  Honert  in  1948  and  extended  by  Cowan  (1965),  namely,  to  treat  the 
flow  of  moisture  from  the  soil  to  the  atmosphere  through  a  plant  as  a 
catenary  process,  i.e.,  limited  by  the  flow  through  the  portion  of  the  path 
with  highest  resistance.  Gardner  (1960)  modeled  the  flow  of  moisture  to  a 
root  with  consideration  of  the  transpiration  demand  rate,  an  approach  extended 
by  Cowan  (1965)  to  include  varying  demand.  These  results  (and  others 
discussed  in  Appendix  B)  established  that  the  concept  of  permanent  wilting 
point  is  questionable  because  plant  wilting  depends  on  both  demand  and 
average  soil  moisture  content. 

35.  An  excellent  discussion  of  three  stages  of  drying  of  a  soil  is 
given  in  Heller  (1968),  while  research  reported  by  Idso  et  al.  (1974)  and 
Jackson  et  al.  (1973)  describes  surface  color  variations  and  soil  moisture 
measurements  related  to  drying  of  an  irrigated  soil.  Theoretical  develop¬ 
ments  relative  to  water  movement  in  drying  soils  have  been  accomplished  by 
Philip  and  de  Vries  (1957)  and  by  Gardner  (1959).  An  interesting  approach 
was  taken  by  Staple  (1974)  who  used  relative  vapor  pressure  of  a  partially 
dry  surface  in  the  Penman  equation  for  potential  evapotranspiration  to  form 
the  upper  boundary  condition  for  flow  of  moisture  within  the  soil. 


soil  characteristic  curves,  and  results  in  an  unbounded  set  of  transition 
curves,  called  scanning  curves,  for  partially  wetted  or  dried  soils  which 
are  then  subjected  to  change  in  the  opposite  direction.  Hysteresis  imposes 
enormous  complexities  for  analytic  solutions  to  problems  of  alternate  wetting 
and  drying.  It  is  handled  in  numerical  solutions  by  using  separate  data 
sets  and  checking  the  progress  of  the  solution  to  see  what  curve  to  use 
for  the  next  step.  New  opportunities  to  deal  with  hysteresis  have  been 
opened  by  Golden  (1980)  using  percolation  theory,  but  like  spatial  varia¬ 
bility,  hysteresis  is  a  problem  yet  to  be  fully  resolved. 

31.  Advances  in  soil  moisture  measurement  techniques  have  been  well 
covered  in  a  review  paper  by  Schmugge,  Jackson,  and  McKim  (1979  and  1980), 
including  discussion  of  both  in  situ  methods  and  remote  sensing  methods  with 
a  listing  of  the  relative  advantages  and  disadvantages  of  each  method. 

Little  would  be  gained  by  a  repetition  of  their  material  here,  but  some 
comments  are  in  order.  In  situ  measurement  is  labor  intensive,  particularly 
when  field  heterogeneity  is  considered,  while  the  remote  sensing  techniques 
are  generally  limited  to  the  surface  layer  of  the  soil,  1  to  5  cm  (1/2  to 
2  in.).  Remote  sensing  of  vegetation  temperature  may  sense  deeper  layers 
since  the  vegetation  temperature  is  dependent  in  part  on  transpirational 
cooling,  but  serious  problems  remain  in  relating  temperature  to  moisture  for 
a  broad  spectrum  of  plant  species  and  environmental  conditions.  One  is  led 
to  the  natural  conclusion  that  remote  sensing  and  in  situ  methods  will  be 
useful  for  soil  moisture  prediction  purposes  only  if  coupled  to  models  which 
more  efficiently  treat  depth  and  spatial  variations  beneath  and  between  the 
measured  soil  volumes. 


Evapotranspiration 


32.  Evaporation  and  transpiration  are  primary  depletion  mechanisms 
for  soil  moisture  and  must  be  considered  in  any  soil  moisture  prediction 
scheme.  Evapotranspiration  is  basically  a  limited  process,  where  the  evapo¬ 
transpiration  rate  is  limited  by  the  supply  of  energy  for  vaporization,  by 
ventilation  to  remove  water  vapor,  or  by  the  water  supply.  It  involves 
solar  energy,  atmospheric  radiation,  sensible  heat  from  the  atmosphere, 
near-surface  winds,  conduction  of  heat  from  the  soil,  moisture  movement 
within  the  soil,  plant  transmission  of  water  to  leaves,  root  growth,  and  the 
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characteristic  curve  that  was  adjusted  to  conform  to  limited  data  and  derive 
representative  curves  of  both  matric  potential  and  hydraulic  conductivity 
for  all  moisture  contents  with  sufficient  accuracy.  Even  without  measurement 
specific  to  a  sample,  moderately  accurate  curves  could  be  generated.  This 
approach  has  been  used  by  several  modelers  (see  Appendix  C)  and  is  consistent 
in  concept  with  the  development  of  tentative  average  relations  for  sand, 
silt,  and  clay  as  done  for  the  current  WES  soil  moisture  prediction  model. 
Examples  of  typical  soil  characteristic  curves  for  USDA  system  soil  types 
are  given  in  Hanks  and  Ashcroft  (1980)  after  work  reported  in  Taylor  and 
Ashcroft  (1972).  Unfortunately,  none  are  known  for  USCS  soil  types,  except 
by  translation  from  USDA  system  soil  types  as  in  Meyer  and  Knight  (1961) , 
although  typical  USCS  soil  curves  could  be  generated. 

28.  A  number  of  analytic  solutions  to  relatively  simple  problems  of 
soil  moisture  prediction  have  been  developed  in  the  past  three  decades.  The 
majority  of  the  solutions  apply  to  problems  of  infiltration  into  dry  soil, 
as  hysteresis  does  not  then  complicate  the  physics  of  water  movement.  While 
of  limited  value  for  extensive  soil  moisture  prediction  in  the  field,  these 
analytic  solutions  are  of  great  value  in  testing  numerical  models,  as  the 
correspondence  of  analytic  solution  and  numerical  solution  should  be  high 
for  identical  system  assumptions. 

29.  Field  heterogeneity  of  soil  properties  has  been  studied  with  some 
success  in  the  treatment  of  spatical  variability.  Field  heterogeneity  is  a 
problem  of  massive  proportions  for  military  traff icability  and  hydrology 
applications,  as  well  as  agricultural  applications,  because  it  reduces  the 
spatial  scale  which  must  be  considered  to  the  order  of  metres,  and  thereby 
increases  the  computational  requirements  by  orders  of  magnitude  relative  to 
consideration  of  the  spatial  scale  of  soil  types.  The  principal  success  to 
date  has  been  in  the  treatment  of  soils  which  have  similar  physical  character¬ 
istics;  i.e.,  they  differ  only  in  scale  at  the  microscopic  level,  but  much 
remains  to  be  done.  Measured  matric  potentials  and/or  hydraulic  conduc¬ 
tivities  have  proven  useful  for  the  estimation  of  the  microscopic  scaling 
parameter,  while  research  on  its  spatial  variation  is  continuing. 

30.  In  addition  to  the  spatial  variation  of  soil  properties,  there 
are  temporal  variations  in  the  relations  between  soil  moisture  content  and 
both  matric  potential  and  soil  moisture  diffusivity  due  to  the  phenomenon  of 
hysteresis.  Hysteresis  requires  the  determination  of  both  wetting  and  drying 


to  military  problems  and  recommendations  for  improvement  of  the  present 
military  approaches  are  presented  in  following  parts.  Additional  details  on 
developments  in  soil  physics,  evapotranspiration,  and  numerical  modeling  are 
given  in  Appendices  A,  B,  and  C,  along  with  a  more  ordered  review  of  the 
literature  of  the  last  three  decades. 


Soil  Physics 

25.  Hydraulic  conductivity  is  a  critical  parameter  of  the  soil  with 
respect  to  soil  moisture  movement  and  its  prediction.  Water  flux  density 
is  given  by  the  product  of  hydraulic  conductivity  and  the  gradient  of 
hydraulic  potential  (Darcy’s  law),  and  enters  all  models  which  deal  with  the 
process  physics  of  soil  moisture  changes.  Hydraulic  conductivity  is  a 
difficult  factor  to  measure  directly;  thus  research  leading  to  methods  for 
calculating  it  from  the  desorption  soil  characteristic  curve  (matric  potential 
versus  soil  moisture  content)  of  a  sample  has  simplified  the  data  acquisition 
burden  for  soil  moisture  prediction. 

26.  The  recent  approaches  to  the  calculation  of  hydraulic  conductivity 
from  a  soil  characteristic  curve  are  revisions  of  the  original  effort  of 
Childs  and  Collis-George  (1950) .  The  method  utilizes  the  soil  characteristic 
curve  as  a  measure  of  the  distribution  of  passages  for  water  movement  within 
the  soil  matrix,  in  contrast  to  prior  efforts  to  model  flow  blockage  due  to 
soil  grains  using  the  grain  size  distribution.  Because  the  soil  characteristic 
curve  is  more  readily  measured  than  hydraulic  conductivity  and  may  be  measured 
using  an  essentially  undisturbed  soil  sample  in  contrast  to  soil  grain  size 
distributions,  the  approach  is  superior  on  two  counts.  Numerous  investigators 
have  extended  the  method,  including  consideration  of  soil  moisture  diffusivity 
(see  citations  in  Appendix  A) .  Results  of  the  method  have  been  found 
sufficiently  accurate. 

27.  Another  development  of  significance  in  reducing  the  data  acquisi¬ 
tion  burden  has  been  the  identification  of  soil  characteristic  curves  for 
typical  soils  representing  the  soil  types  in  the  USDA  system.  Given  knowledge 
of  soil  type  and  a  limited  set  of  site-specific  measurements  of  matric 
potential  versus  soil  moisture  content,  say  at  -15,  -0.06,  and  -0.005 
atmospheres  matric  potential  corresponding  roughly  to  permanent  wilting  point, 
field  capacity,  and  saturation,  respectively,  one  could  use  a  typical  soil 


PART  II:  SOIL  MOISTURE  PREDICTION  IN  AGRICULTURAL  SCIENCE 


22.  Military  or  civilian  applications  of  traf f icability  and  hydrology 
require  knowledge  of  soil  moisture  content  as  a  function  of  depth  for  exten¬ 
sive  areas  in  order  to  estimate  soil  strength  without  in  situ  measurement 
and  to  predict  surface  runoff.  Prediction  of  soil  moisture  content  in  turn 
requires  knowledge  of  precipitation,  surface  flow,  drainage,  evapotranspira- 
tion,  and  soil  type  as  a  function  of  depth  for  the  same  extensive  areas,  and 
for  a  substantial  period  prior  to  the  time  for  which  a  prediction  is  being 
made. 

23.  Severe  limitations  are  imposed  on  soil  moisture  prediction  by  the 
areal  extent  for  which  predictions  are  needed,  the  spatial  and  temporal 
variability  of  precipitation  and  evapotranspiration,  insufficient  data  on 
soil  characteristics  with  depth,  and  inaccessibility  of  some  areas  critical 
to  military  operations.  Even  nominally  homogeneous  areas  of  a  single  soil 
type  have  been  shown  to  have  variations  of  soil  properties  which  must  be 
considered  in  soil  moisture  prediction;  thus  prediction  must  be  made  'or 
fairly  small  areas  or  must  account  for  variability  within  larger  areas. 
Subsurface  soil  type  also  varies  widely  within  nominally  homogeneous  areas 
of  a  specific  surface  soil  type,  and  its  determination  presents  very  serious 
difficulties  due  to  the  large  areas  involved.  A  measure  of  spatial  varia¬ 
bility  of  precipitation  from  cyclonic  storms  and  of  evapotranspiration  may 
be  obtained  from  local  topography,  including  slope  orientation,  while 
evapotranspiration  is  also  influenced  by  vegetation  type.  Convective  storm 
precipitation  is  more  randomly  variable,  however,  and  cannot  be  adequately 
measured  by  surface  networks  of  rain  gages  at  commonly  used  spacings. 

24.  A  coordinated  program  of  measurement  and  computer  modeling  is 
required  to  deal  with  the  problem  of  soil  moisture  prediction.  The  program 
must  be  further  adapted  to  treat  extensive  areas  with  a  useful  degree  of 
accuracy.  Research  in  agricultural  applications  has  led  to  greatly  improved 
understanding  and  definition  of  the  important  factors  for  soil  moisture 
prediction  and  to  vastly  improved  models  and  computational  methods  for  their 
solution  over  the  last  three  decades.  Advances  in  remote  sensing  of  near¬ 
surface  soil  moisture  have  also  occurred.  The  recent  developments  in  agri¬ 
cultural  science  related  to  soil  moisture  prediction  are  discussed  in  this 
part  of  the  report.  A  discussion  of  the  applicability  of  these  developments 


correlations  of  data  in  forming  the  model.  This  historical  point  has  been 
emphasized  because  model  improvement  is  likely  to  require  reconsideration  of 
processes  which  were  lumped  into  the  regression  relationships,  and  because 
subsequent  reports  have  not  communicated  the  early  attempts  to  build  the  model 
from  a  sound  scientific  base,  leading  to  inappropriate  criticism  by  some 
model  reviewers. 

21.  The  WES  soil  moisture  prediction  model  in  current  use  is  therefore 
an  updated  version  of  an  approach  developed  almost  three  decades  ago.  For 
the  time,  the  approach  was  formulated  with  awareness  of  the  known  physics  of 
the  extremely  complex  problem.  For  the  present,  the  model  has  incorporated 
some  increased  knowledge,  but  it  has  been  overtaken  by  models  based  on  more 
complete  analysis  of  problem  details.  Comparison  tests  may  establish  that 
the  current  WES  model  is  still  as  good  as  others  for  traf f icability  and 
mobility  predictions  worldwide.  However,  in  the  opinion  of  the  author,  this 
cannot  last,  and  a  theoretically  based  model  would  meet  with  much  broader 
acceptance. 


actual  moisture  content  of  the  two  layers  just  prior  to  the  precipitation 
event.  Precipitation  in  excess  of  accretion  for  the  two  layers  is  not 
explicitly  treated. 

18.  Seasonally  dependent  depletion  rates  are  used  in  the  model  via 
depletion  relations  for  winter,  summer,  and  transition  (spring  and  fall) 
seasons.  Separate  rates  are  also  required  for  each  soil  type;  sand,  silt, 
and  clay,  and  for  each  layer,  resulting  in  18  depletion  relations  for  the 
model.  The  depletion  relations  for  the  USDA  soil  types;  sand,  silt,  and 
clay,  have  been  correlated  with  soil  types  under  the  Unified  Soil  Classifica¬ 
tion  System  (USCS) ,  also.  The  depletion  relations  were  derived  by  visual 
averaging  of  a  large  number  of  superimposed  graphs  of  moisture  content 
versus  time,  first  to  derive  relations  for  specific  sites,  then  for  several 
sites  with  a  common  soil  type.  The  present  model  uses  polynomials  of  up  to 
sixth  order  that  were  curve-fitted  to  the  averaged  graphs  to  generate 
functions  for  computer  programming  of  the  model.  These  depletion  relations 
are  scaled  to  the  range  of  soil  moisture  contents  between  field  maximum  and 
field  minimum  for  each  application. 

19.  The  model  was  found  to  be  reasonably  accurate  when  tested  against 
field  data  from  over  100  sites  with  fine-grained  soils,  but  as  noted  above 
(paragraph  14) ,  it  was  less  accurate  for  wet  soils  and  soils  with  high 
organic  content  and  clay.  Additional  details  of  the  model  and  comparison  to 
field  data  may  be  found  in  Smith  and  Meyer  (1973)  and  Carlson  and  Horton 
(1959). 

20.  The  early  effort  leading  to  the  present  WES  soil  moisture  model 
was  based  on  knowledge  of  many  atmospheric  and  soil  factors  that  affect  the 
temporal  and  spatial  variation  of  soil  moisture,  particularly  the  work  of 
the  Forest  Service  personnel.  Lull  (1953)  reviewed  the  literature  of  soil 
physics,  plant  physics,  and  evapotranspiration  with  quoted  excerpts  which 
demonstrate  an  awareness  of  the  state  of  knowledge  at  that  time.  Two  progress 
reports  prepared  by  the  Forest  Service  (US  Forest  Service,  1951  and  1952)  also 
include  discussion  of  physical  and  biological  factors  known  to  be  relevant  to 
the  problem  of  predicting  soil  moisture.  Unfortunately,  these  factors  became 
blended  into  brute  force  regression  relations  and  are  no  longer  explicitly 
treated.  Soil  moisture  prediction  proved  to  be  intractable  as  a  problem  in 
physics  at  that  time,  and  it  was  found  necessary  and  sufficint  to  rely  on 


change  from  specific  relationships  derived  for  prediction  sites  to  general 
relations  based  on  soil  type,  this  step  was  necessary  for  model  applications 
to  extensive  areas  without  data.  The  accuracy  retained  indicates  the  basic 
validity  of  the  relations  derived  with  respect  to  dominant  physical  processes, 
but  the  excessive  errors  (judgement  frrtm  Carlson  and  Horton  (1959))  in  the 
prediction  of  soil  moisture  for  poorly  drained  and  wet  sites,  and  for  soils 
high  in  organic  content  and  clay,  indicated  model  limitations.  It  is  also 
notable  that  the  overall  accuracy  level  achieved  was  largely  due  to  good 
predictions  for  low  moisture  contents,  which  are  not  critical  for  traffica- 
bility  or  mobility. 

15.  Additional  studies  in  tropical  areas  met  with  similar  results; 
i.e.,  the  prediction  accuracy  associated  with  USDA  soil  types  rather  than 
with  soil  data  for  particular  sites  was  red-  ced;  however,  the  results  from 
the  model  were  the  best  predictions  available,  and  model  adjustments  for 
improved  predictions  could  be  made  with  the  enlarged  data  base.  After  a 
traf f icability  and  mobiity  program  review  by  a  selected  board  of  consultants 
was  held  at  WES  in  1966  (WES,  1967),  principal  research  emphasis  was  shifted 
to  other  problem  areas.  The  model  was  refined  as  time,  data,  and  support 
allowed,  with  some  proposed  improvements  incorporated  in  the  computer 
program  of  the  model,  written  around  1970  (Smith  and  Meyer,  1973). 

16.  The  WES  soil  moisture  prediction  model  is  a  budget  model  in  which 
a  simple  bookkeeping  procedure  is  used  to  account  for  daily  moisture  changes 
in  two  15-cm  (6-in.)  layers  in  the  upper  30  cm  (1  ft)  of  soil.  (The  model 
has  recently  been  extended  to  90  cm  (3  ft)  with  tentative  relations  subject  to 
validation.)  Soil  moisture  is  increased  by  accretion  due  to  precipitation 
and  reduced  by  depletion  between  precipitation  events.  The  physical  processes 
of  depletion  are  not  explicity  treated.  The  range  of  soil  moisture  variation 
is  held  between  field  maximum  and  field  minimum  moisture  contents  that  have 
been  determined  by  site-specific  measurement  or  by  expressions  derived  from 
multiple  linear  regressions  against  soil  properties  for  similar  soils. 

17.  Soil  moisture  accretion  for  each  layer  is  determined  by  expres¬ 
sions  derived  by  linear  regressions  of  accretion  data  against  precipitation 
or  available  moisture  storage  in  the  layers.  The  former  form  is  used  for 
cases  where  precipitation  is  less  than  available  storage,  while  the  latter 
form  is  used  when  precipitation  exceeds  available  storage.  Available  storage 
is  defined  as  the  difference  between  field  maximum  moisture  content  and 


Soil  property.  A  physical  characteristic  of  the  soil  matrix, 
such  as  porosity,  soil  particle  types,  or  moisture  diffusivity. 

Tensiometer.  An  instrument  for  measuring  matric  potential. 

Trafficability.  The  capacity  of  a  soil  to  withstand  traffic  of 
vehicles. 

Transpiration.  The  loss  of  soil  moisture  or  plant  moisture  by 
vaporization  at  plant  surfaces  after  passage  through  the  plant. 

Tridiagonal  matrix.  A  square  matrix  in  which  the  only  nonzero 
elements  are  on  the  main  diagonal  and  on  the  diagonals  immediately 
above  and  below  the  main  diagonal. 

Weighting  lysimeter.  A  device  for  weighting  a  physically  isolated 
block  of  soil  in  nominally  natural  surroundings  to  measure  evapo- 
transpiration. 

Wetting  curve.  The  soil  characteristic  curve  derived  during 
wetting  of  a  dry  soil  sample. 

Wetting  front.  The  spatial  position  of  the  sharp  gradient  of 
soil  moisture  between  dry  and  saturated  portions  of  a  soil  column 
during  wetting. 

WES  Soil  Moisture  Prediction  Model 

12.  A  brief  review  of  the  WES  soil  moisture  prediction  model  is 
included  in  this  section  to  provide  a  reference  for  readers  who  are  not 
familiar  with  the  model,  and  as  a  guide  to  the  author's  viewpoint  for  those 
who  are. 

13.  Pilot  studies  were  conducted  during  the  period  from  1948  to  1953 
by  WES  and  US  Forest  Service  personnel  to  identify  the  scope  of  the  problem 
of  predicting  soil  moisture  and  to  formulate  a  preliminary  prediction  model 
and  research  plan  for  its  improvement.  Two  preliminary  models  were  developed 
independently  by  the  WES  Trafficability  Section  and  the  US  Forest  Service 
during  this  period.  As  the  Forest  Service  model  afforded  greater  potential 
for  future  development,  it  was  selected  for  continued  effort,  although  each 
model  had  predicted  soil  moisture  about  as  well  for  the  data  available  (WES 
and  US  Forest  Service,  1954,  Vol  3). 

14.  Data  collection  was  extended  via  cooperation  with  other  groups 
and  additional  sampling  by  WES  and  Forest  Service  personnel.  The  enlarged 
data  base  enabled  formulation  of  the  model  for  the  general  soil  types;  sand, 
silt,  and  clay,  from  the  US  Department  of  Agriculture  (USDA)  soil  classifica¬ 
tion  system  (Carlson  and  Horton,  1959).  While  some  accuracy  was  lost  in  the 


51.  The  choice  of  a  particular  set  of  model  characteristics  is  not 
self-evident.  Models  have  been  formulated  with  various  combinations  of 
methods  of  finite  differencing,  nonlinearity  treatment,  hysteresis  treatment, 
and  step  size  selection.  An  optimum  combination  requires  consideration  of 
the  physical  process  treated,  the  intended  use  of  the  results,  resources, 
and  modelers'  points  of  view  concerning  accuracy,  stability,  economy,  and 
convenience.  Direct  comparisons  of  model  performance,  such  as  that  of 
Havercamp  et  al.  (1977)  for  infiltration,  are  required  as  an  objective  guide 
to  model  selection  or  model  development  for  new  applications. 

52.  In  addition  to  models  which  treat  the  movement  of  water  in  the 
soil  matrix  in  one-,  two-,  and  three-dimensional  systems,  several  models 
have  been  developed  which  include  plant  roots  in  the  soil,  and  a  few  include 
explicit  consideration  of  the  entire  plant.  Plant  considerations  are 
important  in  determining  the  depth  to  which  soil  moisture  is  strongly  affected 
by  evapotranspiration,  and  some  redistribution  of  moisture  within  the  soil 
profile  is  effected  by  root  transfer  of  moisture  from  more  moist  to  dryer 
layers. 

53.  Two  basic  representations  of  roots  within  the  soil  volume  are 
used.  A  microscopic  root  model  treats  flow  of  moisture  to  individual  roots 
that  are  considered  straight,  infinitely  long,  and  isolated  from  other 
roots.  The  moisture  uptake  of  the  roots  per  unit  length  is  calculated,  and 
then  multiplied  by  a  measure  of  total  root  length  in  a  volume  to  obtain 
total  water  extraction.  This  representation  of  root-soil  water  interaction 
is  necessary  for  study  of  the  effects  of  moisture  gradients  in  the  vicinity 
of  roots.  A  macroscopic  root  model  treats  flow  of  moisture  to  roots  via  a 
bulk  measure  of  rooting,  such  as  rooting  density.  The  process  physics  of 
flow  to  the  roots  is  not  explicitly  treated  in  these  models,  with  consequent 
loss  of  resolution  and  reduction  of  computational  costs.  Each  approach  has 
been  used  with  success  (see  Appendix  C) . 

54.  Currently  available  soil  moisture  prediction  models,  even  the 
simpler  mass  balance  models  to  be  discussed  in  the  following  section,  are 
limited  by  available  input  data.  Several  papers  have  reported  differences 
between  model  output  and  laboratory  data  which  may  best  be  explained  as  data 
errors,  for  instance,  the  assumption  of  uniform  packing  of  soil  in  a  labora¬ 
tory  test  column.  Moisture  movement  occurs  in  response  to  actual  soil 
packing  density,  while  the  model  must  use  what  the  experimenter  inputs  about 


packing  density.  Field  data  of  sufficient  vertical  and  horizontal  resolution 
for  specification  of  initial  moisture  condition  and  soil  characteristics 
important  to  soil  moisture  prediction  are  extremely  sparse. 

55.  Field  and  laboratory  data  obtained  at  Agricultural  Experiment 
Stations  and  Agricultural  Research  Service  offices  throughout  the  nation  in 
connection  with  various  experiments  and  modeling  efforts  are  frequently 
sufficient  for  modeling  purposes;  however,  no  collection  of  such  data  sets 
was  noted  during  the  literature  review.  A  set  of  desorption  data  for  over 
1800  soils  on  experimental  watersheds  was  published  by  Holtan  et  al. 

(1968,  see  also,  Holtan  et  al.,  1967).  The  data  include  textural  class. 

These  data  were  used  by  Clapp  and  Hornberger  (1978)  to  evaluate  their  power 
curve  representations  of  soil  hydraulic  properties.  The  only  larger  body  of 
data  noted  is  that  gathered  during  the  development  of  the  WES  soil  moisture 
prediction  model  and  reported  with  that  effort. 

56.  Model  data  requirements  vary  greatly,  as  the  application  for 
which  each  model  was  developed  and  the  modeler's  anticipation  of  data 
availability  vary  greatly.  On  the  one  hand,  models  have  been  used  with 
estimated  characteristics  for  typical  soils  and  schematic  (sinusoidal) 
variation  of  evapotranspiration,  while  on  the  other  hand,  some  models 
require  detailed  curves  of  matric  potential  and  hydraulic  conductivity  or 
moisture  diffusivity  versus  soil  moisture  content,  detailed  soil  profile 
data,  and  full  meteorological  data  for  calculation  of  potential  evapotrans¬ 
piration,  as  well  as  precipitation  on  site.  A  separate  research  project 
would  be  required  to  detail  data  required  for  each  model  and  assess  the 
significance  of  each  datum.  While  this  has  not  been  attempted,  inputs  for 
several  of  the  models  reviewed  in  Appendix  C  are  noted  along  with  results  of 
sensitivity  studies  when  reported.  Further  comments  on  data  for  models  are 
included  in  sections  on  applicability  and  recommendations. 

57 .  A  measure  of  cost  for  numerical  modeling  of  soil  moisture  was 
given  by  Wind  and  Van  Doome  (1975).  Their  one-dimensional  model  cost  about 
38c  per  simulation  day  to  run.  It  produced  results  appropriate  to  traffic- 
ability  applications,  including  depth  to  the  water  table,  for  the  top  1  m 

(3  ft)  of  soil.  Programming  convenience  may  result  in  costs  an  order  of 


magnitude  higher  (Havercamp  et  al.,  1977;  Richter,  1980),  while  some  cost 
reduction  may  result  from  tailoring  the  model  to  specific  applications. 
Higher  cost  will  accrue  for  two-  and  three-dimensional  models  with  similar 
vertical  resolution.  Cost  will  also  increase  with  the  number  of  sites  and 
simulation  days  required. 

Mass  Balance  Models 

58.  Mass  balance  (budget)  models  are  based  on  the  law  of  conserva¬ 
tion  of  matter.  They  generally  calculate  soil  moisture  content  of  a  soil 
layer (s)  from  a  prior  moisture  content  by  simple  addition  and  subtraction  of 
moisture  fluxes  due  to  such  processes  as  precipitation,  evapotranspiration, 
runoff,  and  drainage.  Measured,  estimated,  or  otherwise  forecast  precipita¬ 
tion  and  irrigation  are  used  as  inputs,  while  outputs  are  generally  deter¬ 
mined  by  statistical  methods,  including  linear  regressions  based  on  field 
data.  Since  process  physics  is  not  influenced  by  the  character  of  a  model, 
the  seemingly  simple  expressions  used  for  moisture  fluxes  in  mass  balance 
models  must  still  incorporate  physical  processes  via  parameterization  or  in 
the  regression  coefficients.  Darcy's  law  is  not  used  in  these  models,  thus 
rendering  them  tractable  without  numerical  integration  schemes,  but  many  of 
the  models  do  explicitly  treat  evapotranspiration  and  other  aspects  of  the 
total  process  physics  with  separate  expressions.  Most  are  adapted  to 
computer  applications. 

59.  Relatively  simple  one-layer  budget  models  were  developed  by 
Thornthwaite  and  associates  (Thomthwaite  and  Mather,  1954)  and  by  Jensen, 
Robb,  and  Franzoy  (1970).  The  Thornthwaite  model  was  used  for  estimating 
soil  moisture  for  prediction  of  traction  capability.  Moisture  input  was 
given  by  measured  precipitation,  while  both  gravitational  drainage  (above 
field  capacity)  and  evapotranspiration  were  considered  for  moisture  losses. 
Potential  evapotranspiration  via  the  Thornthwaite  method  was  adjusted  for 
soil  moisture  content  below  field  capacity.  The  Jensen  et  al.  model  was 
developed  for  irrigation  scheduling.  It  uses  measured  and  forecast 
precipitation  plus  irrigation  in  soil  moisture  prediction.  Potential 
evapotranspiration  via  the  Penman  method  is  adjusted  by  means  of  crop 
coefficients  to  determine  actual  evapotranspiration  (the  only  loss)  through 
the  growing  season. 


60.  The  two-layer  model  developed  at  WES  has  been  discussed  in  the 
Introduction  of  this  report.  It  is  a  mass  balance  model  with  extensive  use 
of  data  correlation  in  the  model  parameterization. 

61.  Multilayer  mass  balance  models  have  been  developed  by  Baier  and 
Robertson  (1965,  1966),  Jones  and  Verma  (1971),  and  Stuff  and  Dale  (1978). 

Several  additional  multilayer  models  related  primarily  to  crop  yield 
modeling  are  also  noted  by  Hildreth  (1978)  with  a  listing  of  model  equations. 

Baier  and  Robertson  use  their  model  for  latent  evaporation  to  estimate 
potential  evapotranspiration  and  then  combine  potential  evapotranspiration 
with  available  soil  moisture  in  each  layer  to  calculate  actual  evapotrans¬ 
piration.  Their  calculation  includes  an  empirical  coefficient  to  account 
for  soil  and  plant  characteristics  of  each  layer.  Precipitation,  runoff, 
and  percolation  through  soil  layers  are  considered  in  expressions  used  to 
calculate  soil  moisture  change  in  each  layer  each  simulation  day.  Variable 
layer  thickness  is  allowed  in  the  model. 

62.  Hildreth  (1978)  compared  calculations  by  the  Baier  and  Robertson 
model  with  field  data  for  one  month  under  a  wheat  crop.  He  found  the  model 
results  to  be  within  10  percent  of  measured  values  in  six  30-cm  (1-ft)  layers 
at  the  end  of  that  period,  with  the  total  moisture  content  calculated  within 
2  percent  of  the  measured  amount.  He  also  found,  however,  that  several 
combinations  of  model  parameters  resulted  in  similarly  accurate  results,  but 
one  set  of  parameters  resulted  in  poor  model  performance  in  an  independent 
test.  Hildreth's  sensitivity  tests  of  the  model  indicated  no  amplification 
of  errors  in  input  values,  but  output  error  varied  one-to-one  with  input 
error  of  10  percent  in  soil  moisture  capacity,  available  soil  moisture, 
rooting  coefficient,  and  a  soil  dryness  coefficient.  Further,  he  found 
errors  to  be  additive,  leading  to  potentially  large  error  in  model 
calculations. 

63.  The  model  of  Stuff  and  Dale  (1978)  was  developed  to  treat  soil 
moisture  prediction  with  high  water  tables  under  corn.  It  was  based  on 
empirical  relations  developed  from  two  years'  field  data,  including  measured 
water  table  depth  and  capillary  rise  from  the  water  table.  The  model  had 

seven  layers  each  15  cm  (6  in.)  thick.  Evapotranspiration  was  based  on  I 
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measured  pan  evaporation  with  adjustment  factors  for  crop  development  and 
moisture  stress  factors.  Further  model  development  would  be  required  for 
application  to  other  areas  and  to  extend  the  range  of  soil  moisture  used  for 
parameterization,  but  this  model  is  one  of  the  few  mass  balance  models  with 
consideration  of  water  table  effects. 

64.  Bouma  et  al.  (1980a,  1980b)  also  considered  the  effects  of  water 
table  levels  using  several  levels  of  detail  in  their  model  input  data  that 
were  based  on  soil  surveys  and  field  sampling.  They  found  strong  influence 
on  predicted  soil  moisture  from  subsoil  characteristics  and  also  found  that 
using  a  single  field  boring  to  determine  subsoil  type  for  a  given  homo¬ 
geneous  area  on  a  soil  map  was  detrimental  to  predicting  accuracy.  Either 
field  estimates  of  subsoil  type  made  during  soil  surveys  or  multiple  borings 
yielded  superior  predictions.  Their  work  dramatizes  the  importance  of 
subsoil  data,  as  even  the  relatively  simple  prediction  model  they  used  was 
limited  in  prediction  accuracy  by  input  soil  data. 

65.  Schmugge,  Jackson,  and  McKim  (1980),  Kanemasu  et  al.  (1980),  and 
Hildreth  (1978)  should  be  consulted  for  additional  details  and  references  to 
individual  models  of  the  mass  balance  type.  Further  information  is  also 
provided  in  Appendix  C. 


PART  III:  DISCUSSION 


66.  Numerous  advances  in  the  areas  of  soil  physics,  evapotranspira- 
tion,  and  modeling  of  soil  moisture  have  been  made  in  agricultural  science 
since  the  currently  applied  model  for  military  soil  moisture  predictions  was 
developed.  During  that  period  there  have  also  been  several  advances  in  the 
state  of  the  WES  soil  moisture  model  as  additional  data  and  research  have 
been  used  to  improve  the  regression  relations  used;  to  improve  estimation  of 
field  maximum  and  field  minimum  moisture  contents;  and  to  treat  problems  of 
spatial  variability  of  soil  properties,  water  table  influences,  and  problem 
soils  of  high  organic  content.  The  WES  model  has  been  adapted  for  computer 
applications  with  some  revision  of  the  accretion  and  depletion  relations, 
fitted  polynomial  curves  for  these  relations,  and,  in  addition  to  soil  type, 
consideration  of  site  conditions  (by  means  of  a  wetness  index  used  to 
describe  the  influence  of  water  table  and  drainage  conditions  on  potential 
maximum  soil  moisture  content) .  The  model  has  also  been  tentatively  applied 
to  depths  of  1  m  (3  ft)  for  C5A  transport  aircraft  taxiways. 

67.  The  question  of  applicability  of  the  developments  in  agricultural 
science  to  military  problems  of  traf f icability  and  hydrology  is  therefore 
not  simple,  because  the  current  system  is  useful,  tuned  over  many  years  to 
specific  military  considerations,  and  in  place.  In  one  sense,  perhaps  the 
dominant  sense  in  view  of  limited  financial  and  manpower  resources  and 
competing  opportunities,  the  applicability  of  developments  in  agricultural 
science  requires  determination  that  their  implementation  would  result  in 
improvement  to  the  present  system  worth  the  effort.  In  a  less  pragmatic 
sense,  applicability  judgement  would  be  based  on  potential  usefulness  of  a 
specific  development  without  evaluation  of  competitive  stature  with  respect 
to  current  military  methods.  The  following  discussion  strives  for  modera¬ 
tion  and  addresses  opportunities  in  preference  to  certainties. 

68.  In  keeping  with  recent  workshops  on  soil  moisture  prediction 
(Heilman  et  al.,  1978)  and  military  hydrology  (WES,  1979)  it  is  appropriate 
to  consider  short-  and  long-term  research  needs  and  to  consider  amendments 
to  the  present  model  versus  development  of  a  new  approach. 


Short-term  Research  Applications 


WES  soil  moisture  prediction  model 

69.  The  WES  soil  moisture  prediction  model  uses  soil  classification, 
field  maximum  and  field  minimum  soil  moisture  contents,  accretion  and  deple¬ 
tion  relations,  seasonal  transition  dates,  initial  conditions  of  soil 
moisture,  and  precipitation  data. 

70.  Soil  classification.  The  present  model  uses  only  three  soil 
types,  sand,  silt,  and  clay,  combining  several  separately  recognized  soil 
types  in  both  the  USDA  and  USCS  classification  systems  into  these  three 
catagories.  Revision  of  field  moisture  content  factors  and  accretion  and 
depletion  relations  to  incorporate  additional  classes  could  be  done  with  the 
guidance  of  typical  soil  characteristic  curves,  moisture  capacities,  model 
outputs  for  typical  soils,  and  the  guidance  of  process  physics.  Advances  in 
soil  physics  and  modeling  have  provided  bases  for  such  revision. 

71.  Typical  soils  have  been  identified  and  described  only  for  USDA 
types,  while  the  USCS  system  has  been  emphasized  in  military  applications. 

The  above-noted  possibility,  however,  in  conjunction  with  the  association  of 
USCS  soils  with  USDA  soil  types  on  the  USDA  textural  triangle  in  Meyer  and 
Knight  (1961),  may  be  productive  in  the  definition  of  model  parameters  for 
USCS-specif ic  soil  types.  Characteristics  of  typical  soils  and  model 
simulations  to  test  the  proposed  characteristics  of  USCS  typical  soils 
against  field  data  are  direct  applications  of  recent  developments  to 
improvement  of  the  WES  model. 

72.  Research  on  field  heterogeneity  of  soil  properties  has  established 
that  variations  within  nominally  homogeneous  areas  can  be  sufficiently  large 
to  critically  affect  predicted  soil  moisture  content.  On  the  one  hand,  this 
variation  creates  great  difficulty  for  soil  moisure  prediction  models 
because  the  input  soil  characteristics  become  uncertain.  On  the  other  hand, 
the  variability  of  soil  properties  places  bounds  on  the  precision  with  which 
soil  properties  may  defensibly  be  specified  in  field  applications  to  large 
areas.  Judicious  application  of  the  results  of  field  heterogeneity  research 
to  military  problems  may  well  simplify  the  modeling  of  soil  moisture  by 
justifying  use  of  typical  soils  and  moderate  accuracy  models.  The  results 

of  modeling  would,  of  necessity,  be  recognized  as  estimates  with  bounds  of 
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probable  error  stated  as  part  of  the  simulation  output.  In  either  case, 
limitations  on  potential  accuracy  or  justification  of  model  simplification, 
field  heterogeneity  is  applicable  to  military  problems. 

73.  The  modeling  work  of  Bouma  and  associates  is  applicable  to  this 
issue,  in  that  they  have  tested  a  relatively  simple  soil  moisture  prediction 
model  against  several  levels  of  input  data  detail,  and  have  related  this 
work  to  soil  types  defined  by  soil  survey  maps.  Military  applications  of  a 
soil  moisture  prediction  model  to  large  areas  will  necessarily  involve  use 
of  soil  survey  data,  and  both  horizontal  and  vertical  heterogeneity  will 
require  treatment.  This  research  direction  may  lead  to  use  of  different 
soil  types  for  the  two  layers  of  the  WES  model. 

74.  Field  moisture  contents.  Field  maximum  moisture  content  can  now 
be  simply  calculated  for  a  specified  soil  profile  using  one  of  many  numerical 
models  of  infiltration.  The  definition  of  field  maximum  moisture  content 
may  be  used  to  specify  boundary  conditions,  including  water  input.  Summa¬ 
tion  of  soil  moisture  content  in  the  appropriate  soil  volume  after  steady- 
state  infiltration  has  been  reached  gives  the  desired  value.  With  little 
effort  an  existing  model  may  be  adapted  for  this  calculation,  or  values  for 
typical  soils  may  be  precalculated  for  use  in  specific  applications.  The 
cost  of  such  an  approach  would  be  small,  while  these  models  are  so  accurate 
that  significant  differences  between  measured  and  calculated  values  would 
clearly  indicate  that  the  field  situation  was  not  well  known. 

75.  Field  minimum  moisture  content  is  not  as  simple  because  the 
influence  of  plants  and  the  rate  of  potential  evapotranspiration  will  affect 
the  field  minimum  moisture  content  for  a  specified  soil.  Research  in 
evaporation  and  transpiration  may  be  used  to  define  the  gradient  of  moisture 
content  which  would  prevail  for  specified  boundary  conditions  at  the  surface 
and  deeper  in  the  soil  than  the  volume  of  interest,  but  the  potential  for 
significant  improvement  in  modeling  of  subsequent  high  moisture  contents  in 
the  soil  that  are  critical  to  military  operations  is  small.  This  statement 
is  not  true  in  cases  of  high  water  tables  or  significant  subsurface  inflow 
to  a  region  of  interest.  In  such  cases,  subsurface  soil  and  flow  data  are 
required  for  use  in  two-  and  three-dimensional  numerical  models.  See  the 
discussion  of  long-term  research  for  further  comments. 


76.  Accretion  and  depletion  relations.  The  accretion  and  depletion 
relations  used  in  the  WES  model  embody  the  model  process  physics.  These  rela¬ 
tions  incorporate  the  processes  of  infiltration,  redistribution,  drainage,  and 
evapotranspiration.  Since  correlation  relations  such  as  these  frequently 
contain  implicit  error  compensation  in  the  regression  coefficients,  there  is 

a  high  probability  that  accurate,  explicit  introduction  of  one  of  the  processes, 
say  evapotranspiration,  will  not  lead  to  improved  results,  even  after  new 
regression  relationships  are  derived  for  the  residuals.  Short-term  research 
effort  would  be  better  spent  in  other  areas,  such  as  identifying  and 
researching  soil  and  subsurface  flow  conditions  which  violate  the  assumption 
of  a  well-drained  homogeneous  soil  that  is  implicit  in  the  present  model. 

77.  Seasonal  transition  dates.  If  the  possibility  of  frozen  soil  is 
set  aside,  transition  dates  are  probably  dominated  by  changes  in  potential 
evapotranspiration  from  season  to  season.  Several  models  for  calculation  of 
potential  evapotranspiration  are  available,  and  additional  steps  to  convert 
potential  to  actual  evapotranspiration  have  also  been  worked  out.  Short¬ 
term  research  with  a  high  potential  for  success  would  involve  correlations 
of  calculated  potential  and  actual  evapotranspiration  with  accepted  transi¬ 
tion  dates  for  cases  in  the  WES  data  files.  Strong  correlation  is  antici¬ 
pated  between  WES  model  transition  dates  and  periods  of  rapid  transition  in 
calculated  potential  evapotranspiration. 

78.  Initial  soil  moisture  conditions  and  precipitation.  When  measure¬ 
ment  is  not  possible  or  cost  effective,  initial  conditions  may  readily  be 
calculated  by  starting  one  to  three  months  prior  to  the  required  initial 
time,  using  nominal  initial  conditions,  and  running  the  model  for  that 
period.  WES  in-house  research  has  been  done  on  this  problem,  and  many 
agricultural  scientists  have  performed  this  Lest  on  their  models.  It  would 
simply  remain  to  apply  these  previous  results  to  identify  an  optimum  lead 
time  for  given  soil  and  precipitation  conditions. 

79.  Precipitation  timing  and  magnitude,  particularly  in  forecasting 
soil  moisture  contents,  is  a  significant  unknown.  Even  instrumented  areas 
are  seldom  covered  well  enough  to  define  the  spatial  pattern  of  precipita¬ 
tion  from  each  event.  Recent  advances  in  weather  radar  and  the  use  of 

lidar  systems  for  rainfall  measurement  for  calibration  is  frequently  required. 


I 

i 

Because  most  of  the  modeling  work  and  other  research  related  to  soil  moisture 
prediction  in  agricultural  science  has  been  related  to  rainfall  deficiency 
and  irrigation  requirements,  no  major  advances  toward  a  solution  to  the  j 

problem  of  precipitation  were  noted  in  this  literature  review.  < 

New  modeling  approaches  | 

80.  Several  mass  balance  models  have  been  published  which  could  be  ] 
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tested  against  the  current  WES  model  to  determine  their  potential  for  i 

improved  calculation  of  soil  moisture.  Several  of  the  models  are  multilayer 
models  and  may  be  better  able  to  handle  organic  silts  and  moist  areas  than 
the  present  model.  Possibilities  include  replacement  of  the  present  model, 
model  selection  for  different  soil  moisture  conditions,  and/or  development 
of  a  hybrid  model  with  features  of  two  or  more  models  included.  Implementa¬ 
tion  of  these  models  would  be  relatively  easy  due  to  their  simplicity, 
especially  if  supported  by  the  developers*  experience.  Testing  them  to  the 
degree  that  the  WES  model  has  been  tested  would  be  a  larger  effort. 

81.  Preliminary  investigations  of  numerical  models,  particularly 
models  written  in  CSMP  or  another  simulation  language,  would  also  fall  under 
the  category  of  short-term  research.  These  models  work  and  are  directly 
applicable  to  soil  moisture  prediction  problems  for  traf f icability  and 
hydrology,  and  their  utility  for  prediction  of  conditions  which  are  poorly 
handled  by  the  WES  model  needs  to  be  tested.  This  investigation  would  also 
permit  a  critical  evaluation  of  the  WES  model  against  the  state  of  the 

art.  Tests  would  be  most  uieaningful  if  limited  to  cases  critical  to 
traf f icability  or  hydrology,  as  accuracy  at  soil  moisture  contents  so  low 
that  mobility  is  not  impaired  and  runoff  does  not  occur  is  generally  of 
minimum  value  for  military  applications. 

82.  These  newer  models  that  could  be  implemented  with  relatively 
little  investment  may  also  prove  valuable  in  a  two-stage  modeling  approach. 

The  first  stage  would  utilize  the  current  model  to  screen-specific  applica¬ 
tions  for  potential  problems,  e.g.,  rating  cone  indices  between  15  and  125, 
where  it  is  assured  that  impossibility  of  movement  or  certainty  of  movement 
have  been  established  at  these  extremes.  Cases  between  these  limits  would 
then  be  treated  with  a  more  complex  model  to  improve  the  precision  with  which 
potential  mobility  could  be  predicted.  Numerous  models  have  been  developed 
that  are  applicable  to  this  approach.  Tests  are  required  to  determine  whether 
the  approach  would  be  a  cost-effective  improvement  to  the  present  system. 
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Long-term  Research  Applications 

WES  soil  moisture  prediction  model 

83.  This  section  discusses  the  applicability  of  recent  developments 
in  agricultural  science  to  evolutionary  changes  in  the  present  modeling 
approach,  which  is  characterized  by  precipitation-dependent  accretion 
relations  and  moisture  content-dependent  depletion  relations  in  a  mass 
balance  approach.  The  following  section  discusses  revolutionary  changes  in 
approach  to  soil  moisture  modeling. 

84.  The  problems  due  to  field  heterogeneity  of  soil  properties  need 
to  be  addressed  with  field  experiments  and  theoretical  developments.  The 
short-term  approach  outlined  above  may  be  extended  to  incorporate  into  area¬ 
wide  modeling  both  developments  to  the  present  and  future  research  advances 
in  this  area.  Serious  definition  of  the  implications  of  field  heterogeneity 
and  bounds  on  prediction  must  be  addressed  for  multiple-vehicle  operations 
as  well  as  along  a  single-vehicle  path.  Adaption  of  the  current  approach  of 
prediction  of  both  soil  moisture  and  soil  strength  to  area-wide  predictions 
of  mobility  for  straight  and  curved  vehicle  paths  will  require  extensive 
research  just  to  exhaust  present  knowledge  of  the  field  heterogeneity 
problem. 

83.  Adaption  of  the  model  to  enable  effective  use  of  present  and 
developing  capabilities  in  remote  sensing  of  soil  moisture  requires 
consideration  of  a  thinner  surface  layer.  The  physics  of  radiant  emission 
as  a  function  of  water  content  of  the  soil  prevent  sensing  the  moisture 
content  below  the  first  1  to  5  cm  (0.5  to  2  in.).  The  present  top  layer 
thickness  of  15  cm  (6  in.)  in  the  WES  model  is  too  great  to  allow  reliable 
application  of  remotely  sensed  soil  moisture  contents  for  initial  conditions 
or  modeling  accuracy  checks.  Further,  model  predictions  for  heavier  vehicles 
will  require  additional  layers  beneath  the  30-cm  (12-in.)  layer  for  which 
the  model  was  developed,  and  to  which  the  present  data  base  applies  for  the 
most  part. 

86.  Acquisition  of  sufficient  data  to  increase  the  near  surface  depth 
resolution  and  overall  depth  of  the  present  model  would  entail  significant 
expense.  It  would  be  far  more  cost-effective  to  use  numerical  models  to 
generate  simulated  data  for  development  of  accretion  and  depletion  relations. 
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While  this  appears  unwise  on  the  surface,  and  is  likely  to  meet  with  resis¬ 
tance  unless  part  of  a  well  planned  program  which  includes  experimental  data, 
it  is  a  fundamentally  sound  approach  for  several  reasons.  First,  many 
numerical  models  are  very  accurate,  frequently  more  accurate  than  field 
conditions  are  known  for  specification  of  initial  and  boundary  conditions. 
Since  assumptions  as  to  field  conditions  are  involved  in  the  WES  model,  they 
may  be  used  for  initial  and  boundary  condition  specification  for  the  numerical 
model.  Second,  numerical  experiments  are  far  less  expensive  than  field 
experiments  per  simulated  data  point  (versus  "real"  data  point) .  Initial 
and  boundary  conditions  for  the  model  may  be  varied  to  simulate  the  scatter 
of  field  data  and  to  broaden  the  simulated  results  prior  to  performing 
correlations  that  would  produce  simulated  accretion  and  depletion  relations 
for  additional  soil  layers.  Third,  thin  surface  layers  may  be  modeled  to 
correspond  to  remote-sensing  limitations  in  depth,  and  several  numerical 
model  layers  may  be  combined  to  simulate  the  15-cm  (6- in.)  layers  of  the  WES 
model  or  the  limited  resolution  of  neutron  sampling  methods  of  soil  measure¬ 
ment.  Direct  comparison  of  simulation  results  with  field  data  and  WES  model 
predictions  may  therefore  be  made  to  ensure  reliability  of  the  data  set  used 
to  increase  the  vertical  resolutions  of  the  WES  model.  Fourth,  anomalous 
and  unknown  field  conditions  do  not  enter  the  process,  thus  avoiding  con¬ 
tamination  of  the  correlations  with  data  due  to  high  water  tables,  vertical 
variation  of  soil  properties,  or  subsurface  flow.  Fifth,  the  present 
extensive  data  set  upon  which  the  present  model  was  based,  which  contains 
soil  strength  measurements  in  addition  to  soil  moisture  measurement,  could 
be  used  in  expanding  the  depth  resolution  of  the  model. 

87.  The  present  author's  preference  would  be  to  directly  implement  a 
numerical  modeling  approach  to  soil  moisture  prediction.  However,  it  is 
recognized  that  the  present  model  has  several  benefits  related  to  execution 
time  and  simplicity,  as  well  as  being  part  of  the  present  mobility  model. 
Limits  to  defensible  resolution  and  accuracy  and  the  advantages  of  evolu¬ 
tionary  change  mediate  in  favor  of  a  higher  resolution  model  using  the 
present  approach,  but  developed  with  the  assistance  of  the  state  of  the  art 
in  soil  moisture  prediction. 

88.  The  WES  model  does  not  do  well  in  cases  of  organic  soils,  high 
water  tables,  or  subsurface  water  movement.  Specific  field  cases  for  which 


WES  model  performance  has  proven  poor  could  be  examined  via  simulation 
modeling.  Model  inputs  could  be  changed  to  test  hypotheses  as  to  the  cause 
of  WES  model  error  and/or  inaccurate  simulation  by  a  more  complex  model. 
Existing  numerical  models  are  sufficient  to  this  task,  and  improved  defini¬ 
tion  of  the  WES  model  in  response  to  identified  error  sources  would  be  the 
result . 

89.  Topographic  influences  on  soil  moisture  content  are  strong.  Not 
only  is  overland  flow  influenced  by  slope  aspect  and  orientation,  but  also 
evapotranspiration  is  strongly  affected  by  slope  due  to  the  importance  of 
solar  energy  supply  of  latent  heat  for  vaporization.  Studies  have  shown 
that  subsurface  movement  of  water  frequently  leads  to  nonuniform  moisture 
distributions  which  are  not  directly  related  to  surface  relief,  also. 
Potential  evapotranspiration  estimates  using  energy  incident  on  a  sloping 
surface  and  models  of  infiltration  which  consider  runoff  in  two  dimensions 


may  be  applied  to  this  problem,  but  unknown  subsurface  soil  layering  and 
flow  channeling  pose  a  more  difficult  problem  to  solve,  as  specific  applica¬ 
tions  require  extensive  specific  site  data.  Long-term  research  in  the 
similar  problems  created  by  field  heterogeneity  and  topography  is  indicated. 
New  modeling  approaches 

90.  Given  a  decision  to  implement  a  new  modeling  approach  for  either 
traf f icability  or  hydrology  applications  of  soil  moisture  prediction,  a 
large  selection  of  both  mass  balance  and  process-oriented  models  are 
available.  Many  of  these  models  may  be  directly  applied  in  their  present 
form  to  both  watershed  and  soil  strength  applications,  while  long-term 
research  would  doubtless  lead  to  specific  model  adaption(s)  for  optimum  use 
in  strategic  and  tactical  military  activities.  A  two-level  approach  which 
involved  initial  screening  by  an  economical,  fast  model  to  be  followed  by 
more  thorough  modeling  of  situations  that  could  pose  problems  for  military 
operations  is  feasbile  with  present  capabilities.  Superior  performance 
would  be  possible  after  optimization  research  based  on  prior  model 
development . 

91.  The  most  difficult  aspect  of  military  applications  of  soil 
moisture  prediction  is  the  magnitude  of  the  area  involved.  Implementation 
of  models  which  can  use  remotely  sensed  data  (thin  layer  sampling  and  poor 
spatial  resolution)  and  typical  soil  characteristics  appears  necessary. 


Standard  model  treatment  of  variability  of  soil  properties  in  all  three 
dimensions  in  a  way  that  produces  usable  predictions  is  also  necessary. 
Research  results  have  been  published  in  each  of  these  areas,  but  a  concerted 
long-term  research  effort  will  be  needed  to  bring  these  results  to  bear  on 
military  problems  in  an  efficient,  effective  model.  Research  has  also  been 
conducted  on  soil  moisture  prediction  based  on  soil  survey  data.  This 
research  direction  is  also  necessary  to  deal  with  the  area  problem,  since  it 
bears  directly  on  specification  of  typical  soils  for  model  use,  and 
represents  the  most  probable  approach  to  military  field  applications  of  the 
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to  build  state-of-the-art  expertise  at  WES  during  model 
development . 

j>.  Revise  the  current  WES  model  as  defensible,  using  methods 
discussed  in  Part  III  of  this  report,  to  increase  soil 
moisture  prediction  capabilities  during  new  model  development. 
This  recommendation  would  be  particularly  germane,  if  the 
basic  approach  of  the  WES  model  is  selected  for  the  screening 
submodel. 

h.  It  is  strongly  recommended  that  the  model  development  include 
generating  model  capability  to  exploit  remote  sensing  for 
initial  condition  specification,  operational  boundary  condi¬ 
tion  adjustment,  and  simulation  validation.  This  will 
require  treatment  of  a  thin  surface  layer  and  spatial 
averaging  to  simulate  low  spatial  resolution. 

1.  It  is  also  strongly  recommended  that  model  development 
emphasize  the  use  of  typical  soils,  including  a  standard 
procedure  that  involves  supplemental  calculations  for 
alternate  submedium  soil  types  that  are  likely  to  occur 
with  a  given  surface  soil  type. 

X-  In  recognition  of  field  heterogeneity,  uncertainty  of  sub¬ 
medium  soil  type,  unknown  subsurface  flow,  and  precipitation 
variability,  it  is  recommended  that  some  form  of  probability 
forecasting  be  used,  either  for  the  predicted  soil  moisture 
content  or  in  the  form  of  error  estimates  on  the  prediction. 

94.  Numerous  details  have  been  left  out  of  the  steps  relating  to  new 
model  development.  Several  possible  steps  and  details  of  approaches  have 
been  discussed  at  several  places  in  this  report.  It  is  recognized  that  the 
results  of  a  serious  test  of  the  WES  model  against  other  models,  on  cases 
where  a  model  is  likely  to  be  used  and  the  results  are  important,  would  be 
the  best  guide  to  further  effort. 


PART  IV:  RECOMMENDATIONS 


92.  It  if  recommended  that  WES  undertake  an  applied  research  project 
to  develop  an  improved  model  tor  soil  moisture  prediction  for  use  in  military 
hydrology  and  traf f icability  problems  as  soon  as  resources  permit.  It  is 
further  recommended  that  the  soil  moisture  prediction  model  consist  of  two 
basic  submodels  which  are  compatible  but  individually  applicable  to  specific 
problem  areas.  One  submodel  is  envisioned  as  a  screening  model  that  effi¬ 
ciently  identifies  soil  moisture  contents  that  may  lead  to  excessive  runoff 
or  may  impair  mobility  in  military  operations.  The  second  submodel  is 
envisioned  as  a  more  complex,  but  more  accurate,  numerical  model  to  be  used 
for  problem  cases  identified  by  the  screening  model. 

93.  The  following  steps  are  recommended  as  one  orderly  approach  to  an 
improved  soil  moisture  prediction  model  that  includes  staged  development  and 
optimization  based  on  earlier  results: 

ci.  Create  a  computer  file  of  the  extensive  soil  moisture/soil 
strength  data  acquired  in  previous  projects  for  use  in  model 
development . 

b.  Supplement  WES  data  with  data  obtained  from  agricultural 

experiment  stations.  Agricultural  Research  Service  offices, 
and  civil  engineering  departments  engaged  in  soil  moisture/ 
soil  strength  research. 

£.  Select  a  limited  set  of  cases  for  which  field  data  are 

available  that  represent  potential  applications  of  a  military 
soil  moisture  prediction  model. 

d.  Test  the  WES  model  on  the  cases  selected  in  £  using  the  full 
input  data  set  and  data  subsets  representing  probable  data 
limits  in  real  military  operations. 

e.  Contract  representative  soil  moisture  modeling  groups  to  test 
their  models  with  the  same  data  set  and  subsets,  and  to 
report  and  discuss  their  modeling  results  in  a  conference 
format,  including  comparisons  of  accuracy,  cost,  and  model 
versatility  against  the  WES  model. 

This  sequence  of  steps  will  provide  a  factual  basis  for  decision  on  whether 
to  proceed  with  a  serious  attempt  to  modify  or  replace  the  current  model. 

It  will  also  provide  a  basis  for  assessment  of  the  modeling  approaches  most 
likely  to  be  optimum  for  screening  and  detailed  simulation  submodels,  if 
continuation  is  indicated.  Assuming  further  model  development  is  supported 
by  the  facts  collected  in  the  above  steps: 
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APPENDIX  A:  SOME  RECENT  DEVELOPMENTS  IN  SOIL  PHYSICS 


Introduction 


1.  A  number  of  developments  in  soil  physics  which  have  been  cited  in 
connection  with  numerical  modeling  of  soil  moisture  are  briefly  discussed  in 
this  appendix.  Other  developments  which  are  considered  particularly  germane 
to  problems  of  traf f icability  and  soil  moisture  modeling  for  large  and/or 
diverse  areas  are  also  included,  such  as  field  heterogeneity  of  soil 
properties  and  methods  of  calculating  soil  characteristics  curves.  The 
latter  developments  would  be  of  value,  either  to  improve  the  basis  and 
procedures  of  the  current  WES  soil  moisture  model,  or  to  create  a  new  model. 
Recent  development  of  soil  moisture  simulation  models  are  discussed  in 
separate  appendices. 


Fundamentals 


2.  The  fundamental  equations  describing  the  movement  of  water  in 
soils  had  been  derived  by  1952,  about  the  same  time  as  the  basic  concepts  of 
the  WES  model  were  being  formulated  and  implemented.  Effective  solutions  of 
the  equations  required  another  couple  decades  for  additional  analytical  work 
for  special  cases  (for  instance,  by  Philip  and  by  Parlange)  and  for 
numerical  techniques  to  be  developed  (see  Appendix  C) . 

3.  Soil  water  exists  in  gaseous,  liquid,  and  solid  states,  depending 
on  temperature  and  water  content.  The  frozen  case  has  not  been  considered 
in  this  report.  Water  vapor  movement  is  important  during  the  final  stages 
of  soil  drying  and  during  periods  of  strong  temperature  gradients  in 
unsaturated  soils  (Hanks  and  Ashcroft,  1980;*  Taylor  and  Ashcroft,  1972; 
Nielsen  et  al.,  1972;  Philip,  1969;  Hadas,  1968;  Cary,  1966;  Philip  and 

De  Vries,  1957;  and  many  others  cited  in  these).  Temperature  gradients  also 
are  noted  to  affect  liquid  flow  to  a  lesser  extent  than  vapor  flow.  Neither 
factor  is  considered  further  herein  except  as  it  has  been  included  in  a 
model  discussed. 


4.  Movement  of  soil  moisture  as  a  liquid  has  been  described  by  means 
of  Darcy’s  law  and  the  equation  of  continuity  (see  Klute,  1952b  or  Philip, 

1969  for  good  general  discussions,  or  any  soil  physics  text).  Darcy  (1856) 
published  an  empirical  relationship  in  which  the  rate  of  steady  state  water 
movement  in  saturated  soils  was  proportional  to  the  gradient  of  soil  water 
potential,  or 

q  =  KV$  (Al) 

where: 

q  =  the  one-dimensional  flow  velocity  [cm  sec  (volume  per  unit 
time  per  unit  area) 

K  =  the  factor  of  proportionality  known  as  coefficient  of  permeability 

or  hydraulic  conductivity  [cm  sec  ^] 

V  =  the  gradient  operator  [grad  ()=V()=i  - h  j  +  K 

—  ox  — 3y  — 9z 

for  Cartesian  coordinates] 

4>  =  the  hydraulic  potential  [cm]  forcing  the  flow. 

If  hydraulic  potentials  due  to  osmotic  pressures  or  solute  concentrations 
and  gradients  are  ignored,  the  hydraulic  potential,  $  ,  may  be  written  as 
the  sum  of  matric  potential  (soil  water  tension,  or  suction),  V  ,  and 
gravitational  potential,  z,  where  z  is  the  vertical  distance  from  a 
prescribed  datum  level  (arbitrary  for  system)  with  units  [cm]  for  each. 

5.  The  equation  of  continuity  for  this  system  may  be  written  simply, 

|f  =  "V*q  (A2) 

3  -3 

where  0  is  the  volumetric  water  content  [cm  cm  ]  (to  correspond  with  WES 
soil  moisture  model  usage,  note  that  inches  of  water  per  six  inches  of  soil 
corresponds  to  60) .  In  Equation  A2,  convergence  into  a  volume  increases  the 
water  content  with  time,  while  divergence  of  the  water  flow  reduces  water 
content.  Equation  Al  is  now  substituted  for  q  in  Equation  A2  to  produce 

If  =  V.(KVd>)  =  V.(KVf)  +  |f  (A3) 

or  rewritten  for  one-dimensional  flow  in  the  vertical. 


(A4) 


where  $  =  'V  +  z  has  been  substituted  and  simplified  in  the  second  equality. 
Equations  A3  and  A4  are  quite  general;  they  apply  to  homogeneous  or  hetero¬ 
geneous  soils  and  pose  no  special  requirements  on  relations  among  0  ,  $  , 
and  K. 

6.  It  is  convenient  to  introduce  the  concept  of  diffusivity  (Childs 
and  Collis-George,  1950;  and  Klute,  1952a  and  b)  in  order  to  eliminate 
matric  potential  from  Equations  A3  and  A4.  Thus,  defining  the  diffusivity, 

D  ,  as 

D  -  K  ||  (A5) 

and  introducing  D  into  Equations  A3  and  A4  produces,  respectively, 


90 

at 


V • (D  76) 


90_ 
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(A6) 


9t  9z  9z  )  90  *  9z 


(A7) 


where  the  chain  rule  has  been  used  in  each  equation.  It  is  necessary  for 
this  transformation  that  the  soil  be  homogeneous  over  the  domain  of  applica¬ 
tion  of  the  equation,  and  that  K  and  7  be  single-valued  functions  of  0  . 
D  is  then  also  a  single-valued  function  of  0  .  Another  transformation  of 
value  in  many  problems  may  be  made  when  K  and  0  are  single-valued 
functions  of  7  .  It  may  then  be  written  that 


90 

ay 


li-  7.(KVT) 

9 1  V  7  97 


(A8) 


90  91  =  3_  /K  3J\  ,  3K  dj_ 

9 V  9t  ~  9z  9z}  9Y  9z 


(A9) 


While  Equations  A6  and  A7  are  generally  more  tractable  than  Equations  A8 
and  A9,  the  latter  may  be  used  with  ponded  water  on  the  surface,  and  are 
used  in  some  of  the  models  discussed  in  Appendix  C.  According  to  Philip 
(1969)  and  others,  iv.  hards  (1931)  first  set  down  formal  equivalents  of 
Equations  A3  and  A8,  while  Equation  A6  was  first  formulated  by  Klute  (1952a 
and  b)  . 
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7.  The  cautious  reader  will  be  concerned  by  such  extensive  develop¬ 
ment  based  on  an  empirical  relation,  Darcy's  law,  by  the  application  of  this 
law  to  unsaturated  soils  (it  was  derived  from  saturated  flow  data) ,  and  by 
the  requirement  for  single-valued  relationships  among  K  ,  6  ,  and  ¥  . 

Each  concern  is  treated  more  fully  in  the  following  paragraphs. 

8.  The  basic  saturated  form  of  Darcy's  law  has  been  deduced  analyti¬ 
cally  by  Hubbert  (1940)  for  conditions  of  negligible  inertial  forces  in  the 
flow  (low  Reynold's  number  in  soil  pores),  and  Fancher,  Lewis,  and  Barnes 
(1933)  present  data  for  sands  which  indicate  that  the  law  holds  at  Reynold's 
numbers  less  than  unity.  Poulovassilis  (1977)  has  shown  that  even  lower  soil 
pore  Reynold's  numbers  may  be  necessary  to  maintain  Darcian  flow  at  low  water 
contents.  Departures  from  Darcy's  law  have  also  been  observed  for  nearly 
impermeable  clay  materials  (Swartzendruber ,  1968),  and  due  to  inertial 
effects  during  some  transient  flows  (Philip,  1969).  Such  details  are  not 
generally  important  for  soil  moisture  modeling  for  traf f icability  or  military 
hydrology,  however.  For  the  present,  and  for  the  present  purpose,  Darcy's 
law  may  be  considered  valid. 

9.  It  was  assumed  by  Richards  (1931)  during  his  derivation  of  Equa¬ 
tions  A3  and  A8  that  within  the  context  of  soil-moisture  models  Darcy's  law 
holds  for  unsaturated  media.  This  was  confirmed  by  Childs  and  Collis-George 
(1950),  and  by  several  others  since  then.  This  finding  may  be  applied  to 
swelling  soils,  also,  when  flow  is  related  to  the  local  soil  particles  and 
continuity  of  the  total  soil-water  matrix  is  considered  (Philip,  1969). 

10.  The  requirement  for  single-valued  functions  among  volumetric  water 
content,  0  ,  matric  potential,  V  ,  and  water  conductivity,  K  ,  is  not 
normally  met  in  general  soil  water  flows.  It  is  well  met  in  several 
specific  cases,  however,  such  as  infiltration  of  water  into  a  profile 
initially  at  equilibrium  with  a  water  table  or  with  uniform  initial  water 
content,  and  rainage  of  a  saturated  soil.  It  is  also  met  during  periods  of 
wetting  or  drying  of  a  local  soil  volume,  even  when  other  processes  are 
occurring  elsewhere  in  the  soil  column.  In  these  cases,  the  requirements  are 
met  locally,  and  computations  merely  require  consideration  of  whether  the 
soil  has  been  wetting,  drying,  or  unchanging.  The  requirement  is  more 
restrictive  for  analytic  solutions  of  the  equations  than  for  numerical 
solutions;  as  in  the  latter  case,  each  computation  for  steps  in  space  or 
time  is  relatively  independent. 
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11.  The  principal  developments  since  1952  have  been  in  methods  of 
ermining  the  essential  K-6  and  V-Q  relations  for  use  of  the  above 
ations,  formulation  of  analytical  solutions  for  special  cases,  and  in 
dies  of  heterogeneity  of  soil  properties  in  nominally  uniform  areas, 
h  category  is  covered  more  fully  in  the  following  paragraphs. 


Methods  for  Estimation  of  Soil  Hydraulic  Properties 


12.  Childs  and  Collis-George  (1950)  opened  the  way  to  an  entirely  new 
of  opportunities  in  soil  moisture  physics  by  proposing  a  method  for 

culation  of  hydraulic  conductivity,  K  ,  from  the  pore  size  distribution 
a  soil.  Previously,  hydraulic  conductivity  (or  permeability)  had  been 
culated  from  the  particle  size  distribution  via  Kozeny's  (1927)  equation, 
by  equations  derived  by  others  (noted  by  Childs  and  Collis-George) .  The 
r  view  was  that  of  a  distribution  of  passages  for  flow  (pores)  rather  than 
listribution  of  blockages  to  flow  (particles) .  The  subsequent  development 
Millington  and  Quirk  (1961),  Brutsaert  (1967),  Green  and  Corey  (1971), 
ipp  and  Hornberger  (1978),  and  Ahuja,  Green,  and  Chong  (1980)  is 
lommended  reading,  along  with  the  original  paper.  While  each  of  the 
:hod  development  papers  includes  data  for  comparison  to  calculations, 
rticular  emphasis  on  measurement  and  test  of  the  method  was  reported  by 
:kson,  Reginato,  and  Van  Bavel  (1965);  Kunze,  Uehara,  and  Graham  (1968); 
ice  (1972) ;  and  Elzeftawy  and  Mansell  (1975) .  Additional  contributions 
re  been  made  by  Marshall  (1958),  Millington  and  Quirk  (1959),  Jackson  and 
Lsler  (1970),  Rogowski  (1971,  1972a  and  b) ,  Jackson  (1972),  and  Parkes  and 
:ers  (1980).  Emphasis  on  soil  moisture  diffusivity,  D  ,  was  discussed  by 
ite  (1952b),  Bruce  and  Klute  (1956),  and  Poulovassilis  (1977),  for  both 
.culation  of  soil  moisture  movement  and  determination  of  soil  moisture 
iracteristics . 

13.  The  method  uses  the  curve  of  matric  potential,  Y  (soil  water 
ision,  or  suction),  against  volumetric  moisture  content,  0  ,  as  a  measure 
soil  porosity  at  each  water  content.  This  approach  is  much  easier  than 
rempting  to  measure  porosity  directly,  and  possesses  the  advantage  of 

Lng  water-pore  interactions  in  the  primary  data,  t  versus  0  .  The  Y 
:sus  9  data  are  often  taken  only  at  selected  values  of  V  ,  but  typical 


ht  of  site-specific  data,  in  particular  the  stages  of  crop  growth  (crop 
fficient),  the  methods  of  determining  radiant  energy  fluxes,  and  the 
tilation  factor.  However,  while  it  is  clearly  not  a  total  answer  to  the 
blem,  potential  evapotranspiration  has  succeeded  in  bringing  greater 
ceptual  order  to  solutions. 

10.  It  is  interesting  to  quote  from  a  final  report  on  a  project 
itled  "Estimating  Soil  Tractionability  from  Climatic  Data"  (Thornthwaite 

Mather,  1954,  p  401)  in  the  context  of  this  report: 

Similar  computations  for  other  places  and  other  years  support 
the  conclusion  that  soil  moisture  can  be  determined  with  all  needed 
precision  from  climatological  data.  It  is  apparent  from  the  agree¬ 
ment  found  between  measured  and  computed  values  that  the  climatologic 
approach  will  permit  the  accurate  determination  of  the  movement  of 
water  through  soils  and  the  amount  of  storage  in  any  selected  layer 
in  the  soil.  .  .  . 

tory  has  shown  this  to  be  too  strong  a  statement. 

11.  Another  general  approach  to  the  problem  of  evapotranspiration  was 
eloped  by  Baier  and  Robertson  (1965) ,  who  used  single  and  multiple 
ressions  of  evaporation  measured  by  evaporimeter  against  several  atmos- 
ric  inputs  to  derive  an  expression  for  latent  evaporation.  This  approach 

examined  in  the  early  work  on  the  WES  model,  with  similar  correlation 
fficients  (Carlson  and  Horton,  1959  (Appendix  E);  Dortignac  and  Lull, 

1) .  While  correlations  of  moderate  accuracy  may  be  derived,  they  must  be 
pled  to  a  soil  moisture  model  for  prediction.  This  approach  was  not 
lowed  for  the  WES  model,  as  it  requires  detailed  modeling  of  soil 
sture  movement.  The  depletion  relations  developed  were  deemed  accurate 
ugh,  but  were  not  readily  adapted  to  incorporate  evaporation  as  a 
arate  depletion  mechanism.  Further,  the  correlation  coefficients  were 
high  enough  to  inspire  confidence  in  the  approach. 

12.  The  Baier  and  Robertson  model  has  been  related  to  potential 
potranspiration  by  Baier  (1971).  This  model  has  been  used  for  both  soil 
sture  modeling  (Baier  and  Robertson,  1966;  Chieng,  Broughton,  and  Foroud, 
8)  and  crop  models  (Baier,  1973;  Feyerherm,  1977).  It  was  reviewed  and 
ted  for  sensitivity  by  Hildreth  (1977). 

13.  Another  approach  of  personal  interest  to  the  author,  which  has 
n  shown  to  be  valuable  objectively  as  well,  is  that  of  Lettau,  called 
apotranspiration  climatonomy."  The  approach  (Lettau,  1969;  Lettau  and 


p  =  the  air  density 

e  =  the  ratio  of  molecular  weights  of  water  vapor  and  air 

k  =  the  Von  Karman  constant 

p  =  air  pressure 

u  =  the  windspeed  at  the  level  z 

3  3 

z  -  the  specified  height  above  the  surface 

3 

*  the  surface  roughness  parameter  of  aerodynamic  theory 

lg  periods  of  strong  surface  heating  or  cooling,  the  atmosphere  is  not 
rally  stratified,  and  corrections  must  be  made  for  stability  in  Equa- 
B13. 

7.  To  calculate  PET  it  is  necessary  to  measure  or  estimate  H  , 

p  ,  z  ,  T  ,  and  d  .  The  other  factors  in  the  equations  are 
o  a  a  n 

er  constants  or  dependent  on  T  ,  the  air  temperature  at  z  .  This 

3  3 

ntially  reduces  to  determining  the  long-  and  short-wave  radiant  energy 
es  (other  contributions  to  H  generally  being  small),  wind  speed, 
erature,  and  humidity,  as  p  and  zq  may  be  estimated  with  sufficient 
racy  if  necessary. 

8.  The  Thornthwaite  and  Penman  approaches  have  been  compared  for 
opical  area  by  Brutsaert  (1965).  He  found  that  the  temperature-based 
1  was  insensitive  there,  due  to  a  small  annual  variation  of  monthly 
eratures.  The  Penman  equation  produced  evapotranspiration  values 
ing  from  0.71  to  1.16  times  the  measured  values  (weighing  lysimeter) . 
s  et  al.,  (1973)  compared  evaporation  measured  via  lysimeter  with  that 
icted  by  the  Penman  and  van  Bavel  equations  under  conditions  of 
uently  strong  advection  of  sensible  heat.  They  found  that,  during  a 
e-week  period  when  water  was  not  limiting,  the  Penman  equation  predicted 

times  actual  evapotranspiration,  while  the  Van  Bavel  equation  predicted 
times  actual.  Daily  plots  of  hourl>  values  indicate  that  the  latter 
estimate  is  due  to  excessive  response  to  strong  late-af ternoon  winds,  as 
Van  Bavel  results  follow  the  actual  evapotranspiration  much  better  in  the 
ling  hours. 

9.  The  concept  of  potential  evapotranspiration  is  particularly 
iable  because  it  enables  partial  separation  of  the  atmospheric  contri- 
on  from  the  plant  and  soil  contributions  to  the  limits  on  evapotrans- 
ition.  It  has  generally  been  found  beneficial  to  revise  the  equations  in 
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if  it  is  assumed  that  the  exchange  processes  of  heat  and  water  vapor  are 
similar,  LE  and  A  may  be  written  as 

LE  =  -LB  (e  -  e  )  (B5) 

v  o  a 

md 

A  =  -yLB  (T  -  T  )  (B6) 

vo  a 

rfhere 

B^  =  the  water  vapor  exchange  coefficient 

Y  *  the  psychrometric  constant,  required  to  maintain  proper 
units  in  Equation  B6. 

Defining  A  as  the  slope  of  the  e'  versus  T  curve,  it  follows  that 

(T  -  T  )  =  (e'  -  e’)/A  (B7) 

o  a'  o  a 

approximately.  Equation  B7  is  introduced  into  Equation  B6  with  the  result 

A  =  -  (y/A)LB  (e1  -  e’)  (B8) 

v  o  a 

Since  (e,-e,)=(e,-e)+(e  -  e')  ,  it  follows  that  Equation  B8  may 

OH  OH  H  H 

be  written 


A  =  -(Y/A)LBv(e^  -  e&)  +  (Y/A)LBv(e^  -  ej  (B9) 

and  since  potential  evapotranspiration  occurs  from  a  well  watered  surface, 
it  is  appropriate  to  consider  that  surface  saturated,  thus,  e^  =  e^  ,  and 

A  =  (y/A)L(PET)  +  (y/A)LBvda  (BIO) 

where  d  =  (e1  -  e  )  ,  the  vapor  pressure  deficit  at  the  specified  height 
H  H  H 

above  the  surface  noted  earlier.  When  A  from  Equation  BIO  is  substituted 
into  Equation  B4  with  LE  =  L(PET) 

L(PET)  +  H  +  (y/A)L(PET)  +  (y/A)LB  d  =  0  (Bll) 

v  a 

and  solving  for  L(PET)  produces 

(y/A)H  +  LB  d 

l(pet)  — (-Y/i)  +y  (bi2) 

Penman  used  an  empirical  relationship  for  B^  ,  while  Van  Bavel  used  a 
relationship  based  on  neutral  atmosphere  boundary  layer  theory  given  by 


,  _  pek^  _ ^a _ 

7  P  [ln(z  /z  )]2 
a  o 


(B13) 
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ST*  was  then  adjusted  for  a  specific  month  length  (d)  and  average  day 
angth  in  the  month  (h)  by 

PET  =  PET*(h/12) (d/30)  (B2) 

o  derive  the  potential  evapotranspiration,  PET.  Thornthwaite  certainly 
new  better  than  anyone  that  heat  supply  and  ventilation  are  important  to 
vapotranspiration,  but  T  is  more  available  than  any  other  usable  data, 
nd  this  empirical  expression  worked  reasonably  well. 

5.  The  approach  of  Penman  (1948,  1956,  and  1963)  requires  less 
eadily  available  data  for  the  computation  of  potential  evapotranspiration, 
ut  it  is  more  physically  based,  and  has  been  extended  by  Van  Bavel  (1966) 

o  remove  the  empiricism  used  by  Penman.  The  Van  Bavel  development  will  be 
hown  with  note  taken  of  the  differences. 

6.  The  development  begins  with  an  expression  for  evaporation  based 
n  work  of  Dalton  (1802),  where 

E  =  f (u) (e'  -  e  )  (B3) 

O  a 

rhere 

f(u)  =  a  function  of  windspeed  at  some  height  above  the  surface 

1/2 

(Dalton  used  u  ) 

e'  =  the  saturation  vapor  pressure  of  water  at  the  temperature 
of  the  surface,  T 

o 

e  =  the  vapor  pressure  of  the  air  at  that  height 
a 

rhile  the  saturation  vapor  pressure  at  Tq  may  be  found  in  tables,  deter- 
lination  of  Tq  is  more  difficult.  Penman  and  Van  Bavel  thus  sought 
‘limination  of  e^  from  Equation  B3.  To  do  so,  it  is  first  necessary  to 
.nvoke  the  surface  energy  budget, 

LE  +  H  +  A  =  0  (B4) 

rhere 

LE  *  the  evaporation  rate  multiplied  by  the  latent  heat  of 
vaporization 

H  »  all  other  energy  fluxes  at  the  surface,  such  as  solar  radia¬ 
tion,  reflected  solar  energy,  longwave  radiation  to  and 
from  the  surface,  sensible  heat  exchange  with  the  ground, 
and  energy  used  in  photosynthesis 

A  •  sensible  heat  exchange  with  the  atmosphere 

B3 


radiation  and  heat  exchange  with  air  and  ground)."  The  value  of  the  concept 
was  affirmed  by  Penman  (1948,  1956,  and  1963)  in  his  formulation  of  a 
potential  transpiration  equation.  Nearly  all  models  of  evapotranspiration 
used  in  current  numerical  models  for  soil  water  depletion  derive  from  either 
Thornthwaite' s  or  Penman's  approach,  often  with  modifications,  but  using 
an  estimate  of  potential  evapotranspiration  as  a  basis  for  calculating  an 
actual  transpiration. 

4.  Potential  evapotranspiration  is  a  function  of  meteorological 
factors  with  small  contribution  from  soil  heat  flux,  thus  it  would  appear 
that  this  relationship  could  be  exploited  to  derive  a  means  of  calculating 
the  former  from  the  latter.  Realizing  that  data  for  a  straightforward 
approach  to  the  problem  would  frequently  be  lacking,  Thornthwaite  (1948) 
formulated  his  equations  for  potential  evaportranspiration  in  terms  of 
station  monthly  average  temperature  and  day  length  only.  Penman,  on  the 
other  hand,  chose  to  use  a  combination  method  using  a  Dalton  (1802)  style 
equation  of  evaporation  in  terms  of  ventilation  and  vapor  pressure  gradients, 
and  the  balance  of  energy  fluxes  at  the  evaporating  (or  transpiring) 
surface.  These  models  are  presented  in  greater  detail  in  the  following 
paragraphs.  While  the  input  data  requirements  are  greater  for  the  Penman 
equation,  they  are  not  excessive.  The  Penman  equation  and  numerous  revised 
equations  derived  from  it  are  generally  more  accurate  (Pelton,  King, 
and  Tanner,  1960;  Tanner  and  Pelton,  1960;  Van  Bavel,  1966  -  good  revised 
equation;  Hanks  et  al.,  1973;  and  many  others).  Budyko  (1956)  has  developed 
a  similar  approach.  Thornthwaite  (1948)  formulated  his  theory  empirically  in 
terms  of  station  temperature.  He  defined  potential  evapotranspiration  for  a 
standard  12-hour  daylight  period  and  30-day  month,  PET*  in  units  of  cm  of 
water,  as 

PET*  =  1.6  (10  T/I)m  (Bl) 


A 


where 


T  =  the  monthly  mean  station  air  temperature,  °C 
I  =  a  heat  index  for  the  station  given  by 


I 


i; 


i  =  (T/5) 


1.514 


>  V_  A  V  . 


w 

O 


'A 


% 


B2 


APPENDIX  B:  SOME  RECENT  EVAPOTRANSPIRATION  DEVELOPMENTS 


Introduction 

1.  A  number  of  recent  developments  in  the  theory  and  modeling  of 
evaporation  and  transpiration  are  discussed  in  this  appendix.  Work  was 
selected  for  review  in  large  part  on  the  basis  of  its  use  in  numerical 
modeling  of  soil  moisture.  Theoretical  and  experimental  developments  are 
emphasized  here,  as  treatment  of  evaporatranspiration  in  numerical  models  is 
discussed  in  Appendix  C. 

2.  Symon's  (1867)*  has  said  that  "evaporation  is  the  most  desperate 
branch  of  the  desperate  science  of  meteorology."  While  some  progress  has 
been  made  in  the  ensuing  century,  our  understanding  of  evapotranspiration  is 
still  insufficient.  The  difficulties  encountered  in  the  treatment  of 
evapotranspiration  accrue,  in  part,  from  the  fact  that  it  is  not  simply  a 
meteorological  problem.  There  are  distinct  contributions  to  actual 
evapotranspiration  from  soil  and  plant  physics  as  well.  Basically, 
evapotranspiration  may  be  viewed  as  a  limited  process,  in  the  sense  that  the 
rate  of  water  loss  is  limited  by  available  energy  and/or  the  rate  of  water 
movement  to  the  ground  surface  and/or  the  rate  of  water  movement  through 
plants  and/or  the  rate  of  vapor  removal  by  the  atmosphere.  The  virtual 
intractability  of  the  problem  of  evapotranspiration  derives  from  the 
complexity  of  each  of  these  limiting  processes  and  their  interaction. 

Atmospheric  Limits  to  Evapotranspiration 


3.  Perhaps  the  most  valuable  concept  put  forth  in  the  field  of 
evapotranspiration  is  that  of  potential  evapotranspiration,  set  down  by 
Thornthwaite  in  the  Report  of  the  Committee  on  Transpiration  and  Evapora¬ 
tion,  1943-44  (Thornthwaite  et  al.,  1944).  Potential  evapotranspiration  has 
been  defined  in  several  slightly  different  ways,  but  an  acceptable  defini¬ 
tion  is:  "that  evapotranspiration  rate  which  occurs  from  a  well  watered 
surface,  i.e.,  the  rate  limited  only  by  meteorological  factors  (including 


*  References  are  collected  at  the  end 


by  Warrick,  Mullen,  and  Nielsen  (1977)  with  good  results  in  both  cases.  The 
spatial  variation  of  1  is  presently  under  study  to  find  measures  of  the 
spatial  applicability  about  a  measured  point.  Autocorrelation  functions  and 
other  statistical  techniques  appear  to  be  important  in  this  effort  (personal 
communication,  D.  R.  Nielsen  and  J.  Wagenet,  1981). 

32.  Vertical  heterogeneity  has  been  studied  analytically  in  the  case 
of  evaporation  from  a  water  table  through  nonhomogeneous  soils  by  Hadas  and 
Hillel  (1972)  and  is  discussed  with  a  particular  example  of  stochastic 
heterogeneity  by  Philip  (1980) .  Study  of  infiltration  into  layered  soils 
was  the  object  of  one  of  the  earlier  numerical  models  by  Hanks  and  Bowers 
(1962) ,  and  many  of  the  models  discussed  in  Appendix  C  have  been  applied  to 
the  vertically  heterogeneous  case. 

33.  Other  papers  related  to  spatial  variability  discovered  during  the 
literature  search  are:  Cassel  and  Bauer  (1975)  -  a  study  of  spatial  varia¬ 
bility  below  tillage  depths  in  fields;  Cameron  (1978)  -  a  study  of  varia¬ 
tions  of  the  soil  water  characteristic  curves  and  calculated  hydraulic 

2  2 

conductivity  from  five  sites  at  six  depths  in  a  225  m  (2400  ft  )  plot; 

Bell  et  al.  (1980)  -  a  study  of  spatial  variability  of  surface  moisture  from 
58  large  field  sites  of  400' X  400  m,  16  hectares  each  (1310  X  1310  ft,  40 
acres  each)  for  use  in  evaluation  of  remote  sensing  techniques;  and  two 
studies  of  watersheds  by  Rogowski  (1972b)  and  Peck,  Luxmoore,  and  Stolzy 
(1977).  These  references,  and  the  others  cited  earlier,  should  be  consulted 
for  additional  information  and  work  cited  in  their  references. 

34.  A  corollary  problem,  that  of  predicting  soil  moisture  from  soil 
survey  or  soil  map  data,  is  discussed  in  Appendix  C.  It  is  likely  that 
heterogeneity,  like  hysteresis,  is  a  problem  which  must  be  treated  in  any 
model  which  is  widely  applied,  as  it  also  is  simply  a  fact  of  life  in  the 
soil-water  domain. 


29.  In  addition  to  horizontal  and  vertical  spatial  heterogeneity  of 
soil  properties,  the  prediction  of  soil  moisture  is  also  affected  by  spatial 
variations  of  precipitation,  other  weather  events,  ground-water  movement, 
and  the  temporal  variation  of  the  same  influences.  The  WES  model  requires 
measured  local  precipitation  (as  do  other  models)  for  useful  accuracy,  for 
instance.  Fortunately,  some  progress  is  being  made  on  the  problems. 

30.  Philip  (1980)  describes  two  basic  forms  of  heterogeneity,  deter¬ 
ministic  and  stochastic  (see  also  the  discussion  by  Freeze  (1975)).  Deter¬ 
ministic  heterogeneity  refers  to  cases  where  soil  properties  vary  in  a  known 
way;  whereas,  stochastic  heterogeneity  refers  to  random  variability.  Deter¬ 
ministic  heterogeneity  includes  the  case  of  scale-heterogeneity,  where 
similar  media  are  involved  and  the  variation  may  be  identified  in  variation 
of  a  local  microscopic  length  scale.  According  to  Warrick,  Mullen,  and 
Nielsen  (1977),  similar  media  have  identical  porosities  and  the  same 
relative  particle  and  pore  size  distributions.  Stochastic  heterogeneity 
also  includes  a  simpler  form  in  which  the  statistics  of  variation  are 
independent  of  location  or  time.  The  simpler  forms  are  yielding  to  research 
to  some  extent.  The  following  discussion  focuses  on  developments  for 
similar  media. 

31.  Similarity  theory  for  porous  media  was  introduced  by  Miller  and 
Miller  (1956) .  Experimental  work  by  Klute  and  Wilkinson  (1958)  and 
Wilkinson  and  Klute  (1959)  supported  the  concept  with  data  for  soil  water 
characteristic  curves  (f  versus  0  ),  hydraulic  conductivity  values,  and 
infiltration  flow  into  similar  media.  Philip  (1967,  also  1980)  proposed 
that,  if  two  media  differ  geometrically  only  in  their  characteristic  length 
scales,  say  1^  and  1  ,  then  matric  potential  and  hydraulic  conductivities  of 
the  two  media  may  be  related  by 

1lH'l  =  V2  and  K1/;L12  =  K2/X2  *  (A11) 

The  values  of  1^  for  several  soils  may  be  determined  from  either  measured 
matric  potentials  or  measured  hydraulic  conductivities.  The  method  was 
tested  on  several  soils  ranging  from  a  clay  to  a  fine  sand  by  Reichardt, 
Nielsen,  and  Biggar  (1972a)  and  found  to  hold.  This  texture  range  is 
extreme,  compared  to  the  assumptions,  and  is  encouraging  for  more  general 
application  of  the  approach.  The  method  was  also  tested  by  Russo  and  Bresler 
(1980)  and  applied  to  three  data  sets  containing  several  hundred  observations 


some  approximation,  many  have  achieved  excellent  correspondence  with  field 
or  laboratory  measurements.  Recent  exploration  of  the  applicability  of  a 
branch  of  abstract  probability  theory,  percolation  theory,  to  soil  moisture 
problems  by  Golden  (1980)  has  opened  new  opportunities  to  deal  with 
hysteresis;  however,  much  more  work  is  still  to  be  done. 


Field  Heterogeneity  of  Soil  Properties 


27.  Variation  of  hydraulic  and  physical  soil  properties  throughout  a 
field  of  nominally  identical  soil  type  (field  heterogeneity)  places  often 
severe  restrictions  on  the  accuracy  of  any  prediction  technique  applied  to 
the  entire  area.  Philip  (1980),  in  an  excellent  overview  of  the  subject, 
refers  to  the  "enormity"  of  the  problem,  meaning  not  "enormousness,"  but 
rather  "monstrous  wickedness,  deviation  from  normal  type,  that  which  is 
abnormal,  a  monstrous  offense."  Field  heterogeneity  generates  staggering 
problems  in  the  general  case,  but  several  approaches  to  solution  of  real 
problems  have  been  partially  successful  because  the  basic  physical  character 
of  a  nominally  homogeneous  field  (one  soil  type,  for  instance)  is  similar 
from  point  to  point. 

28.  Spatial  variation  of  soil  properties  has  been  noted  throughout  the 
development  of  the  WES  soil  moisture  prediction  model  (WES,  1951  and  1952; 

WES  and  US  Forest  Service,  1954;  Carlon  and  Horton,  1957  and  1959;  Collins, 
1971;  and  Broadfoot  and  Burke,  1958),  and  a  special  study  was  conducted  in 
the  area  near  WES  during  1958-1960  by  US  Forest  Service  and  WES  personnel  to 
measure  variability  for  guidance  in  model  development  and  application 
(Carlson  and  McDaniel,  1967).  For  four  loess-derived  soil  series,  they 
found  greater  variability  within  series  than  between  series,  and  concluded 
that  minor  topographical  variations  were  important.  Nielsen,  Biggar,  and 
Erh  (1973)  measured  several  soil  hydraulic  and  physical  properties  for  six 
depths  at  twenty  locations  in  a  150-hectare  (370-acre)  field  in  California. 
The  field  had  been  graded  some  time  earlier  for  improved  irrigation  effi¬ 
ciency,  and  had  been  disturbed  in  the  upper  60  cm  (2  ft)  by  farming  opera¬ 
tions.  They  found  variations  of  steady  infiltration  rate  from  0.5  to  45.7  cm 

day  ^  (0.2  to  18.0  in.  day  ^)  and  variation  of  steady  hydraulic  conductivi- 

-1  2  -1  -1 
ties  trom  about  10  to  roughly  10  cm  day  (0.04  to  40  in.  day  ).  Other 

strong  variations  of  properties  are  noted  in  the  reference  cited. 


and  other  porous  media  by  Gardner  (1959)  ;  to  the  relations  of  external 
conditions  to  drying  of  soils  by  Gardner  and  Hillel  (1962) ;  to  bare  soil 
evaporation,  drainage,  and  storage  by  Black,  Gardner,  and  Thurtell  (1969); 
to  redistribution  of  irrigation  water  by  Gardner,  Hillel,  and  Benyamini 
(1970);  and  to  power  series  solutions  of  the  flow  equation  by  Scott  et  al. 
(1962).  Equations  for  calculation  of  matric  potential  from  volumetric  water 
content  and  hydraulic  conductivity  from  water  content  used  by  Campbell  (1974) 
and  Clapp  and  Hornberger  (1978)  (cited  above)  were  based  on  the  Gardner, 
Hillel,  and  Benyamini 's  (1970)  developments. 

24.  One  of  the  powerful  analytical  techniques  in  solving  the  diffusion 
equation  is  to  reduce  the  number  of  variables  by  integration  of  an  assumed 
function  of  one  of  the  variables  which  reasonably  approximates  the  real 
situation.  This  assumed  shape  may  be  iterated  to  provide  more  accurate 
expressions.  This  approach  has  been  taken  by  Parlange  (1972) ,  Aylor  and 
Parlange  (1973) ,  and  Parlange  and  Braddock  (1980)  for  problems  of  one¬ 
dimensional  infiltration,  infiltration  into  layered  soils,  and  an  accurate 
approximate  solution  of  the  diffusion  equation. 

25.  The  method  of  Green  and  Ampt  (1911)  has  been  described  by  Philip 
(1974)  as  a  primitive  integral  method  using  a  step  function  for  the 
advancing  moisture  profile.  As  with  the  Boltzman  transformation,  integral 
methods  reduce  the  partial  differential  equation  to  an  ordinary  differential 
equation,  which  is  easier  to  solve.  The  Green  and  Ampt  approach  has  been 
used  with  varying  modifications  by  Bouwer  (1969)  for  infiltration  into 
nonuniform  soil,  Raats  (1973)  for  the  analysis  of  unstable  wetting  fronts, 
James  and  Larson  (1976)  for  modeling  infiltration  and  distribution  of 
intermittent  water  applications,  and  Jarrett  and  Fritton  (1978)  for  the 
analysis  of  entrapped  air  effects  on  infiltration. 

26.  Various  other  authors'  analytical  solutions  or  techniques  have 
been  used  in  development  of  numerical  models,  including  the  study  of  one¬ 
dimensional  infiltration  into  homogeneous  soil  by  Fok  and  Hansen  (1966) , 
infiltration  and  runoff  for  small  plots  by  Swartzendruber  (1974)  and 
Swartzendruber  and  Hillel  (1975),  vertical  infiltration  by  Bruts^’rt  (1977), 
infiltration  in  layered  soils  by  Takagi  (1960)  and  by  Fok  (1970),  infiltra¬ 
tion  into  crusted  soils  by  Hillel  and  Gardner  (1969)  and  by  Ahuja  (1973), 
and  even  flow  in  deformable  porous  media  by  Narasimhan  and  Witherspoon 
(1977).  While  each  analytical  solution  to  a  "real"  problem  has  involved 


problems  and  exact  solutions  which  may  be  used  to  evaluate  numerical 
techniques.  In  this  regard,  analytical  solutions  are  preferable  to  field 
data,  as  the  required  answer  is  precise.  Both  analytical  and  numerical 
techniques  must  be  checked  against  field  data  in  the  final  analysis,  however. 

21.  The  analytical  approach  has  been  led  by  E.  C.  Childs,  J.  R.  Philip, 
W.  R.  Gardner,  and  more  recently,  J.  Parlange  during  the  last  three  decades. 
Only  work  related  to  the  numerical  modeling  effort  has  been  referenced  in 
this  report. 

22.  The  work  of  Childs  on  the  problem  of  predicting  hydraulic  conduc¬ 
tivity  for  soils  has  been  noted  above.  Philip  (1957a  and  b)  opened  a 
series  of  papers  on  the  theory  of  infiltration,  which  is  particularly  suited 
to  analytical  treatment  since  problems  may  be  formulated  with  uniform 
antecedent  conditions  and  the  soil  is  raontonically  wetted.  Many  problems 

of  interest  may  be  treated  as  horizontal  infiltration,  thus  converting 
Equation  A6  to  the  simpler  diffusion  equation.  Philip  (1969)  has  written  a 
thorough  review  of  the  theory  of  infiltration,  while  Philip  (1975a  and  b) 
has  presented  an  analysis  of  stability  during  infiltration.  He  has  also 
reviewed  progress  in  the  solution  of  nonlinear  diffusion  equations  (Philip, 
1974),  which  is  a  valuable  confirmation  of  the  fact  that  analytical  solu¬ 
tions  may  be  developed  for  important  problems.  Philip  also  contributed  to 
the  early  application  of  numerical  techniques  to  soil  moisture  problems, 
as  noted  in  Appendix  C. 

23.  Gardner  and  colleagues  have  utilized  a  Boltzman  transformation, 

i/2 

y  =  x/ (DQt)  ,  to  reduce  the  horizontal  diffusion  equation,  90/3t  = 

3  (D  •  30/3x)3x  ,  from  a  partial  differential  equation  to  an  ordinary 
differential  equation. 
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with  consequent  simplification  of  the  problem.  They  have  also  written  the 

diffusivity  in  terms  of  the  volumetric  moisture  content  0  as  D  =  Dq  exp 

6(0  -  0  )  on  the  basis  of  reasonable  fit  to  measurements  and  earlier  work 
o 

by  Wagner,  which  also  enables  analytic  solution  for  some  problems.  In  these 
transformations,  x  is  the  horizontal  distance,  Dq  is  the  diffusivity  at 
0q  ,  t  is  time,  and  B  is  a  constant.  This  approach  was  applied  to  one¬ 
dimensional  infiltration  by  Gardner  and  Mayhugh  (1958);  to  dr’  z  of  soils 


f  versus  K  and  0  versus  K  relations  has  been  found  inadequate  for 
glass  bead  media  (Topp  and  Miller,  1966) ,  for  sandy  loam  soil  (Topp,  1969) , 
and  for  fine  sand  (Vachaud  and  Thony,  1971).  Hysteresis  calculation  remains 
important,  nonetheless,  as  it  is  imperative  that  it  be  considered  in  applica¬ 
tion  of  the  equations  of  soil  moisture  movement,  and  it  is  not  readily 
measured.  The  relative  magnitude  of  41  -  0  variations  due  to  hysteresis 

and  field  variability  has  been  studied  by  Royer  and  Vachaud  (1975).  They 
conclude  hysteresis  must  be  considered  in  field-wide  applications,  as 
V  -  0  variations  due  to  hysteresis  are  greater  than  standard  deviations  due 
to  spatial  variability  reported  by  Nielsen,  Biggar,  and  Erh  (1973).  The 
above  authors  who  have  measured  41  versus  K  hysteresis  conclude  that  it 
is  small  enough  to  ignore  in  most  applications. 

18.  Numerous  analytical  and  numerical  models  of  soil  moisture  flow 
have  been  formulated  with  consideration  of  hysteresis.  This  consideration 
increases  the  computation  time  and  input  data  requirements  of  the  models, 
but  hysteresis  is  simply  a  fact  of  life  in  the  soil-water  domain. 


Analytical  Solutions  for  Special  Cases 


19.  Analytical  solutions  to  the  soil  moisture  flow  equations  are  of 
great  value,  but  they  are  rare  because  of  the  difficulty  of  solving  nonlinear 
partial  differential  equations,  such  as  Equations  A3,  A6,  and  A8.  The  value 
of  analytical  solutions  derives  from  the  exactness  of  the  solution  (within 
the  limits  of  the  solution),  the  general  simplicity  of  calculations  for 
specific  times  or  locations,  and  the  potential  of  learning  about  the 
fundamental  structure  of  the  equation  solution.  Numerical  methods 
(discussed  in  Appendix  C)  are  capable  of  treating  more  general  problems,  but 
each  application  is,  in  a  sense,  separate  from  the  others. 

20.  The  complexity  of  most  applications  of  soil  moisture  flow  equa¬ 
tions  for  soil  moisture  prediction  related  to  traf f icability  or  hydrology 
far  outstrip  the  current  analytical  solutions  of  the  equations.  The  need 
for  answers  can  largely  be  met  by  numerical  methods  by  application  of  brute 
force  solutions  via  digital  computers,  and  this  is  likely  to  remain  true  for 
some  time.  The  fundamental  equations  are  simply  too  difficult  to  solve  in 
general  under  the  conditions  of  the  water-soil  system.  The  analytic  solu¬ 
tion  contribution  will  be  greatest  in  the  provision  of  precisely  stated 


were  required.  Such  applications  are  feasible  using  typical  curves  for  soil 
types,  while  increased  accuracy  should  follow  from  shifting  the  curves  to 
correspond  to  site-specific  measurements  while  retaining  their  characteristic 
shapes.  This  approach  is  similar  to  that  developed  using  tentative  average 
depletion  relations  for  sand,  silt,  and  clay  in  the  current  WES  model 
(Carlson  and  Horton,  1959) . 

16.  A  phehomenon  of  importance  in  the  determination  of  typical  curves 
and  use  of  any  soil  hydraulic  property  data  is  hysteresis.  Hysteresis  is  a 
system  property  in  which  the  state  of  the  system  and  changes  of  that  state 
depend  on  the  system  history.  The  relation  between  matric  potential  and 
volumetric  water  content  exhibits  this  property,  which  is  also  well  known 
from  magnetism  studies.  During  the  process  of  drying,  water  is  removed  from 
pores  in  the  soil,  but  the  removal  is  resisted  by  capillary  forces  at  the 
small  neck  adjacent  to  an  air-filled  pore.  When  the  matric  potential 
decreases  below  a  critical  value  equal  (but  opposite  sign)  to  the  capillary 
force,  the  pore  suddenly  drains  completely  because  capillary  forces  are 
lower  in  the  larger  diameter  pore.  During  wetting,  however,  the  pore 
gradually  fills  up  to  its  maximum  cross  section  at  less  negative  matric 
potentials.  When  the  water  reaches  decreasing  cross  sections,  the  pore 
suddenly  fills  completely  due  to  increasing  capillary  forces  for  smaller 
cross  sections.  Thus,  the  pore  may  range  from  full  to  empty  at  a  given 
matric  potential,  depending  on  whether  wetting  or  drying  is  occurring.  This 
point  was  noted  earlier  in  the  discussion  of  single-valued  relations  among 

V  ,  0  ,  and  K.  Any  soil  physics  text  may  be  referenced  for  further 

discussion  and  figures.  Nielsen  et  al.  (1972)  and  Hanks  and  Ashcroft  (1980) 
were  referenced  for  the  above. 

17.  Hysteresis  is  known  to  be  important  in  the  Y  versus  6 
relations,  and  it  has  been  shown  to  exist  in  the  Y  versus  K  relations. 
While  numerous  experimental  determinations  of  hysteresis  in  soil  hydraulic 
properties  have  been  made,  they  are  far  less  common  than  determinations  of 
the  drying  curves  of  4*  versus  0  due  to  experimental  difficulties. 

Studies  of  hysteresis  in  undisturbed  field  samples,  or  in  the  field,  are 
even  less  frequent.  An  attempt  to  formulate  a  model  for  calculation  of 
hysteretic  Y  versus  0  curves  by  Poulovassilis  (1962)  and  basically  used 
by  Poulovassilis  and  Tzimas  (1974  and  1975)  for  study  of  hysteresis  in  the 


curves  of  ¥  versus  0  have  been  determined  for  all  soil  types  in  the  USDA 
system  (none  are  known  for  USCS  soil  types,  except  by  translation  from  the 
USDA  type  data  as  in  Meyer  and  Knight  (1961) .  Examples  are  given  in  Hanks 
and  Ashcroft  (1980)  after  work  reported  in  Taylor  and  Ashcroft  (1972) . 

These  typical  curves  may  be  anchored  for  a  particular  soil  at  the  measured 
values  of  V  versus  0  ,  and  thereby  provide  continuous  data  for  the  soil. 
Moisture  tension  data  at  0.005  atm,  0.06  atm,  and  15  atm  matric  potential 
were  reported  for  the  soils  used  in  development  of  the  WES  soil  moisture 
prediction  model,  and  additional  data  at  0.33  atm  and  3  atm  matric  potential 
are  also  reported  for  some  soils  (Meyer,  1976).  Soil  studies  conducted  at 
Agricultural  Experiment  Stations  located  at  Land  Grant  universities  through¬ 
out  the  nation  will  provide  an  extensive  data  base  for  particular  soils  in 
addition  to  examples  published  in  various  journal  papers.  Empirical  expres¬ 
sions  for  the  41  versus  0  relation  have  also  been  derived  by  Gardner, 
Hillel,  and  Benyamini  (1970) ,  Campbell  (1974) ,  and  Clapp  and  Hornberger 
(1978).  The  latter  paper  reports  comparison  of  soil  hydraulic  properties 
calculated  via  their  expressions  with  data  for  1,800  soils  reported  by 
Holtan  (1967) .  Regression  models  for  soil  hydraulic  characteristics  in 
terms  of  soil  physical  properties  have  been  developed  by  Gupta  and  Larson 
(1979) ;  the  use  of  physical  methods  to  improve  soil  designation  with  respect 
to  drainage  properties  has  been  discussed  by  Bouma  (1973) . 

14.  Several  equations  of  generally  similar  structure,  but  with 
differing  assumptions  and  parameter  values,  have  been  presented  for  the 
calculation  of  hydraulic  conductivity  versus  moisture  content  from  curves  of 
matric  potential  versus  moisture  content  in  the  above-cited  references. 
Brutsaert  (1967),  Green  and  Corey  (1971),  Jackson  (1972),  Rogowski  (1972a), 
and  Gardner  (1974)  each  discuss  the  various  approaches  and  should  be 
consulted  in  lieu  of  further  discussion  herein. 

15.  The  significance  of  these  development  is  perhaps  obvious,  but  it 
will  be  stated  for  emphasis.  The  application  of  the  equations  for  soil 
water  relations  developed  in  the  above  section  requires  knowledge  of  4* 
versus  0  ,  K  versus  0  ,  and/or  K  versus  4*  at  each  point  in  the  soil 
modeled.  Application  of  these  equations  in  analytical  or  numerical  models 
to  the  extensive  --ireas  of  water  sheds  or  areas  of  military  operations 
requiring  mobility  predictions  would  be  impossible  if  site-specific  data 


Baradas,  1973;  K.  Lettau,  1974;  Hall  1977;  Lettau,  Lettau,  and  Molion, 

1979)  uses  a  simplified  measure  of  soil  water  content  and  includes  both 
atmospheric  and  soil-plant  effects  through  parameterization  (see  also 
Appendix  C) .  It  also  utilizes  both  absorbed  global  radiation  (calculated 
via  shortwave  radiation  climatonomy)  (Lettau  and  Lettau,  1969)  and  precipita¬ 
tion  as  joint  forcing  functions  for  evapotranspiration.  The  concept  of 
potential  evapotranspiration  is  therefore  not  necessary,  nor  are  the  many 
adjustments  (such  as  crop  coefficient),  although  this  physically  based  model 
must  incorporate  these  realities  in  the  parameterization.  This  specific 
approach  is  not  yet  widely  accepted. 

14.  An  evapotranspiration  model  developed  by  Ritchie  (1972)  has 
received  relatively  wide  acceptance  and  is  frequently  referenced.  It  is 
based  on  the  Penman  approach,  but  incorporates  revisions  for  treating  a 
growing  row  crop.  Kanemasu,  Stone,  and  Powers  (1976)  tested  the  Ritchie 
model  with  field  data  from  crops  of  soybean  and  sorghum,  finding  that  daily 
and  seasonal  estimates  of  evapotranspiration  agreed  with  lysimetric 
measurements.  Another  relatively  well  used  approach  is  that  of  Shaw 
(1963),  a  model  for  estimation  of  soil  moisture  under  corn.  He  uses  pan 
evaporation  data  for  Iowa  as  a  measure  of  evapotranspiration,  which  is  then 
adjusted  with  a  crop  growth  factor  and  water  stress  factor  for  plant  response 
to  high  evaporation  demand.  These  factors  are  taken  from  figures  in  the 
paper.  Approximately  80  percent  of  the  estimates  of  June  and  August  soil 
moisture  were  within  0.5  in. /ft  over  the  5-ft-deep  profile.  This  model  was 
partially  an  outgrowth  of  earlier  work  (Denmead  and  Shaw,  1959) .  In  later 
work  (Saxton,  Johnson,  and  Shaw,  1974a  and  b)  emphasis  was  shifted  toward 
the  combination  method  of  Penman.  The  Penman  approach  with  modifications 

for  the  specific  problems  of  data  availability  and  objectives  of  the 
application  has  been  used  by  Jensen,  Wrigh'.,  and  Pratt  (1971);  Hanson 
(1976);  and  Morton  (19/5).  Tn  the  latter  paper,  specific  attention  has  been 
paid  to  a  higher  estimate  of  potential  evaporation  that  would  occur  in  a 
large  area  where  water  was  limiting  (thus  reducing  evapotranspiration  and 
increasing  air  temperature  and  heat  advection) . 

15.  The  Penman  equation  has  been  demonstrated  to  be  reasonably  valid 
for  periods  of  a  day  or  hour,  but  the  Thornthwaite  equation  must  not  be  ust  1 
for  periods  much  shorter  than  a  month.  Martin,  Worm,  and  Wilson  (1979) 


noted  the  failure  of  the  Thornthwaite  method  to  follow  weather  events  over 
even  periods  of  one  week  in  their  comparison  of  several  evaporation  pan 
types,  but  inexplicably  conclude  that  it  may  be  used  for  a  daily  moisture 
balance  technique  for  irrigation  scheduling. 

16.  Several  micrometeorologically  oriented  methods  of  estimating 
evapotranspiration  have  been  developed  during  the  last  three  decades  which 
do  not  use  utilize  the  concept  of  potential  evapotranspiration.  They  are 
formulated  to  calculate  evapotranspiration  from  measurements  directly,  using 
the  energy  balance,  the  energy  balance  in  conjunction  with  the  Bowen  ratio 
(sensible  heat  flux  to  atmosphere/latent  heat  flux  [evaporation]),  or 
aerodynamic  properties  of  the  atmospheric  boundary  layer.  These  methods  are 
discussed  by  Tanner  (1967),  as  well  as  several  other  authors.  As  they 
generally  require  much  more  detailed  data  than  the  above  techniques  and 
are  basically  unsuited  to  potential  traf ficability  and  hydrology  models, 
they  will  not  be  discussed  further  herein. 


Soil  Limits  to  Evaporation 


17.  The  atmospheric  properties  control  evapotranspiration  under  many 
conditions  when  water  is  not  limiting,  as  discussed  above.  When  water  is 
limiting,  however,  the  evaporation  rate  from  the  surface  and  transpiration 
rate  from  vegetative  surfaces  is  limited  by  soil  and  plant  physical  (and 
chemical)  properties.  Much  of  the  progress  in  dealing  with  soil  limitations 
to  evapotranspiration  has  come  from  numerical  modeling  of  the  process. 
Several  papers  dealing  with  experiment  or  analysis  without  numerical 
techniques  have  been  selected  for  brief  discussion  here,  as  they  emphasize 
evapotranspiration  limiting  by  soil  water. 

18.  Three  stages  of  the  drying  of  soils  have  been  identified  with 
respect  to  water  availablility  (Idso  et  al.,  1974).  The  first  stage  is  that 
of  unlimited  water  for  evaporation,  the  second  stage  is  a  falling  rate  stage 
where  evaporation  is  limited  by  vapor  diffusion  from  relatively  free  water 
beneath  the  air-dry  surface,  and  the  third  stage  is  limited  by  adsorptive 
forces  holding  water  to  individual  soil  particles.  These  three  stages  have 
been  discussed  qualitatively  for  decades,  including  by  US  Forest  Service  and 
WES  personnel  involved  in  development  of  the  WES  soil  moisture  prediction 
model,  but  they  had  been  demonstrated  experimentally  only  in  laboratory 


studies  before  1974.  The  approach  used  by  Idso  et  al.  included  measurement 
of  the  surface  reflectivity  for  shortwave  radiation  (albedo)  as  it  changed 
with  surface  soil  moisture  content.  The  change  from  the  first  to  the  second 
stage  of  drying  was  particularly  evident  via  this  technique.  Generally,  the 
surface  would  become  air  dry  during  the  afternoon  of  one  day,  but  became 
remoistened  via  upward  transport  of  water  during  the  night.  It  would  then 
become  air  dry  again,  but  at  an  earlier  hour  of  the  day.  Clearly,  the 
division  between  stages  one  and  two  is  dependent  on  the  rate  of  water  loss, 
as  well  as  soil  factors. 

19.  The  above  results  were  obtained  during  a  series  of  intensive  soil 
measurements,  as  reported  by  Jackson  et  al.  (1973).  They  sampled  at  depths 
of  0.5,  1,  2,  3,  4,  5,7,  and  9  cm  (six  replicates)  every  half  hour  for  the 
full  24  hours  of  sixteen  days  during  the  first  38  days  after  field  irriga¬ 
tion.  Data  for  surface  flux  was  obtained  via  a  weighing  lysimeter.  The 
data  for  four  days  are  plotted  in  three-dimensional  perspective  figures, 
which  clearly  demonstrate  the  complexity  of  near-surface  water  fluxes. 

While  many  previous  experiments  had  been  conducted  under  laboratory 
conditions,  very  few  had  been  conducted  to  show  diurnal  variation  in  the 
field.  In  a  following  paper,  Jackson  et  al.  (1974)  compare  their  measured 
soil  water  fluxes  with  values  calculated  by  the  theory  of  Philip  and  De  Vries 
(1957) .  They  found  the  theory  predicts  water  fluxes  best  at  intermediate 
water  contents,  but  some  serious  deviations  between  measurement  and  theory 
occur  at  high  or  very  low  water  contents.  One  point  raised  was  the  need  for 
very  accurate  values  of  the  moisture  diffusivity  for  application  of  the 
theory.  This  paper  brings  out  the  mutual  limits  to  accuracy  created  by 
approximations  required  for  analytical  solutions  and  experimental  diffi¬ 
culties  in  precisely  measuring  important  soil  physical  properties. 

20.  A  discussion  of  the  three  stages  of  drying  and  the  movement  of 
the  bone  dry  front  in  soils  is  presented  by  Heller  (1968)  with  an  excellent 
overview  of  the  physics  of  saturated  and  unsaturated  moisture  flow.  It  is 
recommended  reading,  particularly  for  its  thoroughness  and  conceptual 
clarity.  The  influence  of  surface  residue  and  evaporation  potential 
(provided  by  ventilation  and  infrared  lamps)  was  studied  in  laboratory  soil 
columns  by  Bond  and  Willis  (1970)  with  emphasis  on  the  effect  on  the  first 
stage  of  drying,  while  Hanks,  Gardner,  and  Fairbourn  (1967)  studied  separate 


effects  of  wind  and  radiation  on  evaporation  rate  in  laboratory  columns  via 
a  similar  exposure  to  heat  lamps  and  fans.  Their  results  show  strong  surface 
cooling  under  ventilation  and  surface  heating  under  radiation.  The  total 
evaporation  from  the  surface  was  similar,  with  slightly  greater  totals  under 
ventilation  for  two  of  the  three  soil  types  tested.  The  first  stage  of 
drying  is  also  stressed  in  an  experimental  technique  developed  by  Arkin, 
Ritchie,  and  Adams  (1974)  for  measurement  of  the  effect  of  surface  mulches 
on  evaporation  from  the  surface.  Their  technique  has  further  utility,  in 
that  it  may  be  used  to  measure  evaporation  from  the  bare  surface  of  cropped 
fields  to  provide  data  for  separating  evaporation  from  transpiration  when 
total  evapotranspiration  has  been  measured  or  estimated.  This  separation  is 
important  in  crop  yield  modeling,  and  may  be  important  in  soil  moisture 
modeling  to  separate  deep  depletion  by  roots  from  more  shallow  depletion  via 
surface  evaporation. 

21.  The  control  of  evaporation  rates  by  water  movement  in  the  soil 
has  been  studied  theoretically  and  experimentally  by  W.  R.  Gardner  and 
colleagues  (Gardner,  1959;  Gardner  and  Hillel,  1962;  Black,  Gardner,  and 
Thurtell,  1969).  The  theoretical  work  was  based  on  the  diffusion  form  of 
the  basic  nonlinear  partial  differential  equation  of  moisture  flow  using  an 
exponential  approximation  for  the  soil  water  diffusivity  (see  discussion  in 
Appendix  A) .  They  have  aided  formation  of  a  sound  theoretical  basis  for  the 
above  discussion  of  soil  limited  evaporation.  Soil  limited  evaporation 
enables  one  to  calculate  evaporation  from  soil  moisture  data  under  a 
limited  set  of  conditions,  but  the  approach  cannot  handle  the  case  of  fully 
wet  surfaces  or  shortcircuiting  of  the  soil  moisture  path  by  plant  roots. 

This  approach  was  also  used  by  H.  R.  Gardner  (1973)  to  analyze  experimental 
data  on  evaporation  from  laboratory  columns  of  water  additions  in  several 
amounts  with  differing  evaporation  intervals.  The  total  water  added  was 
constant  in  three  separate  tests,  as  was  the  total  time  for  evaporation. 
Departures  of  calculated  total  evaporation  from  measured  under  three  condi¬ 
tions  of  water  amount  and  evaporation  interval  were  small,  but  errors  of  up 
to  18  percent  were  noted  for  individual  daily  values.  The  effect  of  crust 
formation  by  rain  versus  flooding  applications  of  water  was  studied  by 
Bresler  and  Kemper  (1970)  to  clarify  the  importance  of  surface  crusts  on 
infiltration  and  subsequent  evaporation.  It  was  found  that  the  crust  formed 
by  rain  reduces  total  evaporation  by  approximately  25  percent.  A  combined 


application  of  the  developments  in  the  study  of  soil  limits  on  evaporation 
is  reported  by  Staple  (1974).  He  modified  Penman's  equation  by  including 
the  relative  vapor  pressure  of  a  partially  dry  surface  soil,  then  used  this 
modified  equation  to  provide  the  upper  boundary  condition  for  solving  the 
flow  equation  within  the  soil.  He  cites  several  investigators  who  had 
previously  used  atmospheric  vapor  pressure  as  the  upper  boundary  condition, 
but  his  mating  of  atmosphere  and  soil  water  is  far  better.  He  found  reason¬ 
able  correspondence  between  measurments  from  small  cylinders  of  soil  imbedded 
in  a  fallow  plot  and  calculations,  but  again  problems  were  encountered  with 
both  the  limiting  assumptions  and  data  requirements  of  the  theory,  and 
difficulties  in  field  sampling  of  the  soil. 


Plant  Limits  to  Transpiration 


22.  The  third  component  of  the  evapotranspiration  system  is  the 
plant.  The  plant  exerts  both  passive  and  active  influences  on  transpiration 
in  systems  of  interest  for  traf f icability  or  hydrology.  The  passive 
influence  derives  from  the  plant's  presence  in  both  the  soil  and  the  lower 
atmosphere  through  distributed  root,  branch,  and  leaf  systems.  The  roots 
access  water  from  deeper  soil  layers  than  is  available  for  surface  evapora¬ 
tion,  while  the  branches  and  leaves  enhance  the  exposure  of  evaporating 
surfaces.  The  plant  also  exerts  active  influence  on  the  stream  of  transpira¬ 
tion  water  at  both  the  roots  and  the  leaves.  Additional  supplies  of  soil 
moisture  for  transpiration  are  made  available  by  root  growth,  both  deeper 
.nto  the  soil,  and  also  horizontally  into  untapped  soil  volumes.  Resistance 
to  vapor  diffusion  from  leaf  stomata  is  increased  as  moisture  stress 
increases  due  to  a  change  of  shape  by  stomate  guard  cells,  thereby  reducing 
transpiration.  These  passive  and  active  plant  roles  increase  survival 
potential  for  the  plant  under  varying  moisture  conditions,  but  they  also 
greatly  increase  the  complexity  of  the  processes  of  soil  moisture  depletion. 

23.  It  is  desirable  in  the  applied  disciplines  of  engineering  and 
agriculture  to  simplify  the  problem  of  soil  water  removal  by  plants  by 
setting  upper  and  lower  limits  of  soil  moisture  which  is  available  for 
transpiration.  The  upper  limit,  field  capacity,  is  set  by  gravity  drainage 
of  soil  water,  while  the  lower  limit,  permanent  wilting  point,  is  set  by 
permanent  wilting  of  the  plant  and  the  consequent  cessation  of  transpiration. 
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The  quantity  of  water  available  between  these  two  soil  moisture  contents  is 
defined  as  available,  or  extractable,  water.  With  these  approach,  transpira¬ 
tion  occurs  until  the  available  water  is  used  up,  then  ceases.  Unfortunately, 
the  natural  system  is  not  that  simple,  and  departures  of  the  real  system 
from  the  simplistic  system  cannot  be  ignored  in  most  cases. 

24.  Three  dominant  process  details  that  prevent  use  of  the  simplistic 
field  capacity/permanent  wilting  point  model  of  soil  water  depletion  by 
plants  are: 

a.  Soil  water  moves  under  the  influence  of  moisture  potential 
gradients,  thus  the  moisture  content  at  the  surface  of  a  root 
is  lower  than  in  the  bulk  of  the  soil. 

b.  The  rate  of  water  loss  varies  strongly  in  time  in  response  to 
both  atmospheric  demand  (potential  evapotranspiration)  and 
the  state  of  growth  of  the  plant. 

c.  Gravity  drainage  of  water  from  many  soils  continues  several 
days  beyond  the  2-day  period  normally  used  in  the  definition 
of  field  capacity;  thus  water  which  would  have  drained  from 
the  soil  (or  evaporated)  in  the  absence  of  vegetation  has 
been  partially  utilized  for  transpiration. 

25.  A  valuable  organizing  concept  for  the  study  of  transpiration  was 
put  forth  by  Van  den  Honert  in  1948.  He  treated  the  movement  of  water  for 
transpiration  through  the  plant/soil  system  as  a  catenary,  or  chainlike, 
process.  In  a  catenary  process  the  overall  rate  is  limited  to  the  rate  of 
the  slowest  partial  process;  that  is,  it  is  a  series  system.  This  organizing 
concept  permits  the  study  of  individual  partial  processes  under  the  assump¬ 
tion  that  other  partial  processes  are  not  limiting.  In  this  way  it  is 
similar  to  the  study  of  potential  evapotranspiration  under  the  assumption 
that  water  is  not  limiting  at  the  soil  surface. 

26.  Gardner  (1960)  analyzed  the  flow  of  soil  water  to  individual 
roots  under  different  soil  moisture  conditions  and  transpiration  rates,  with 
emphasis  on  dynamic  effects.  His  analytic  solution  to  the  flow  equation  for 
water  movement  to  a  root  in  cylindrical  coordinates  enabled  calculation  of 
soil  moisture  gradients  in  the  vicinity  of  the  root  for  different  average 
moisture  contents  and  for  different  transpiration  rates.  The  bulk  soil 
moisture  content  at  the  wilting  point  for  the  plant  based  on  the  soil  matric 
potential  at  the  root  surface  varied  significantly  as  the  transpiration  rate 
was  varied.  Thus,  wilting  would  occur  at  higher  bulk  moisture  content  under 
strong  transpiration  demand  when  steep  moisture  gradients  formed  in  the 
vicinity  of  the  root. 


27.  Gardner  and  Ehlig  (1963)  carried  the  analysis  further  with  more 
accurate  assumptions  about  the  plant  response  to  moisture  stress,  based  on 
recent  measurements.  They  found  the  transpiration  rate  to  be  nearly  a 
linear  function  of  soil  moisture  percentage  when  below  a  threshhold  value. 
Their  analysis  led  also  to  the  conclusion  that  soil  water  remains  available 
for  transpiration  at  matric  potentials  well  beyond  the  15-bar  value  normally 
used  as  permanent  wilting  point.  The  limits  to  transpiration  rates  due  to 
soil  water  movement  to  the  roots  under  steady  evaporative  demand  when  water 
movement  occurred  between  soil  layers  and  regions  of  high  versus  low  rooting 
density  were  also  considered. 

28.  Cowan  (1965)  extended  the  prior  work  with  a  more  complete  treat¬ 
ment  of  the  transpiration  stream.  He  noted  first  that  the  catenary  process 
description  put  forth  by  van  den  Honert  did  not  consider  parallel  paths  of 
flow  through  the  system.  The  principal  parallel  path  for  soil  water  transfer 
to  the  atmosphere  consists  of  soil  surface  evaporation  combined  with  the 
through-plant  transpiration  stream.  Even  within  the  latter,  the  multipli¬ 
city  of  roots,  xylem  vessels,  branches,  and  leaves  leads  to  a  complex 
series/parallel  network  through  which  the  flow  occurs.  Cowan  therefore 
formulated  an  analytical  model  of  transpiration  which  included  series/ 
parallel  paths  in  the  plant  environment  and  the  plant  itself. 

29.  Cowan  also  noted  that  diurnal  variations  of  potential  transpira¬ 
tion  strongly  affect  the  plant-soil  moisture  response,  and  that  considera¬ 
tion  of  steady  evaporative  demand  was  not  sufficient  to  analyze  the  entire 
problem  of  transpiration  limits.  He  therefore  introduced  diurnal  variation 
to  his  model  and  analyzed  the  response  of  soil  moisture  content  and  matric 
potential  over  several  diurnal  periods.  His  results  indicate  reduced 
transpiration  during  diurnal  periods  due  to  soil  moisture  gradients  near  the 
roots,  despite  adequate  bulk  soil  moisture  content. 

30.  Denmead  and  Shaw  (1962)  measured  the  transpiration  rates  of  corn 
plants  in  a  field  experiment  to  investigate  the  effects  of  varying  potential 
transpiration  under  a  range  of  soil  moisture  conditions.  They  found  that 
actual  transpiration  fell  below  potential  transpiration  despite  "adequate" 
bulk  soil  moisture  content  under  conditions  of  high  potential  transpiration. 
At  potential  rates  of  6  to  7  mm/day  (0.24  to  0.28  in. /day),  the  actual 
transpiration  began  to  fall  relative  to  potential  transpiration  at  soil 


APPENDIX  C:  NUMERICAL  MODELING  OF  SOIL  MOISTURE 


Introduction 

1.  This  appendix  presents  an  ordered  review  of  the  majority  of  soil 
moisture  modeling  efforts  over  the  past  three  decades.  The  discussion  is 
organized  along  a  combination  of  historical  and  modeling  group  lines,  where 
the  set  of  publications  for  each  coherent  group  of  investigators  is  briefly 
reviewed  in  the  order  of  their  initial  publication.  As  neither  library 
resources  nor  available  time  were  unbounded  in  this  project,  the  following 
review  could  not  be  exhaustive,  but  it  is  more  complete  than  many  others 
dicovered  during  the  course  of  the  effort. 

2.  While  this  appendix  can  provide  an  ordered  guide  to  the  litera¬ 
ture,  an  outline  of  model  characteristics,  and  identification  of  laboratory 
and  field  data  used  for  model  verification,  it  cannot  prove  Lie  superiority 
of  one  model  over  another,  detail  all  models  published,  nor  detail  the 
published  comparisons  of  simulation  results  with  data.  Moreover,  since  the 
characteristics  of  a  best  model  are  contradictory,  e.g.,  convenience  versus 
computational  speed  and  cost,  any  attempt  at  absolute  judgement  would  be 
absurd. 

3.  Each  of  the  soil  moisture  models  discussed  in  this  appendix  has 
been  published  in  a  research  journal  after  peer  review.  As  a  result  of  this 
prior  filtering  process,  each  of  the  approaches  is  "good"  in  some  signifi¬ 
cant  sense,  and  poor  correspondence  with  data  is  seldom  noted.  In  fact,  the 
model  limitations  seldom  result  from  inaccuracy  of  the  numerical  method 
used,  as  precise  analytical  solutions  can  be  reproduced  to  any  accuracy  for 
which  time  and  money  are  available.  The  model  limitations  derive,  rather, 
from  the  completeness  and  accuracy  with  which  all  relevant  processes  have 
been  included  (or  are  known).  Similarly,  validation  of  the  models  is 
complicated  by  the  lack  of  certain  knowledge  of  the  physical  processes 
relevant  to  particular  measurements.  In  many  respects  the  models  more 
accurately  reflect  their  input  assumptions  and  data  than  do  the  results  of 
laboratory  and  field  experiments.  Departures  of  simulated  results  from 
real  data  are  therefore  often  indicators  of  omission,  and  are  often  clues  by 
which  improved  understanding  through  additional  study  is  provoked. 


4.  An  excellent  discussion  of  numerical  modeling  appropriate  to  the 
papers  reviewed  in  this  appendix  is  contained  in  the  book  by  Remson, 
Hornberger,  and  Molz  (1971).  Hillel  (1977)  has  compiled  a  summary  of  modeling 
effort  which  indicates  the  range  of  possibilities  for  model  applications, 
although  with  nearly  exclusive  emphasis  on  the  approach  he  uses.  Excellent 
comparisons  of  analytic  and  numerical  models  as  applied  to  one-dimensional 
infiltration  were  conducted  and  reported  by  Havercamp  et  al.  (1977)  and  by 
Vauclin  et  al.  (1977).  In  these  papers,  several  different  basic  approaches 
to  the  problem  were  used  and  their  simulations  were  compared  to  an  analytic 
solution  and  laboratory  data.  Freeze  (1969)  also  reviewed  numerical  models 
for  soil  moisture  flow  which  had  been  published.  This  paper  contains  a 
large  table  of  model  characteristics,  enabling  ready  comparison  of  the 
models.  Each  of  these  works  is  recommended  reading. 


5.  Several  terms  are  used  in  the  following  model  discussions  which 
have  been  defined  in  the  main  text  of  this  report  and  discussed  in  Appendix  A. 
A  general  discussion  of  their  significance  is  presented  here,  as  an  introduc¬ 
tion  to  the  model  reviews. 

6.  Richards  (1931)  combined  the  equation  of  continuity  of  water  sub¬ 
stance  with  Darcy's  law  to  derive  the  Richards  equation  used  by  most  of  the 
model  developers  as  a  starting  point.  Darcy's  law  states  that  moisture  flux 
is  proportional  to  the  gradient  of  hydraulic  potential,  while  the  equation 
of  continuity  simply  relates  the  change  of  water  content  of  a  volume  due  to 
net  flow  across  the  volume  boundaries.  The  Richards  equation  is  written  in 
one  dimension  (vertical)  as 


39_  =  3_  ,  3Y  3K 
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where 

9  =  volumetric  moisture  content 
t  =  independent  variable  of  time 
z  =  independent  variable  of  depth 

K  =  hydraulic  conductivity,  the  coefficient  of  proportionality  in 
Darcy's  law 

T  =  matric  potential,  with  hy  .raulic  potential  H  =  T  +  z 
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7.  This  equation  poses  several  inconveniences  for  direct  solution, 
i  the  specific  approaches  used  by  various  modelers  to  deal  with  the 
jation  lead  to  several  different  models.  First,  the  equation  is  written 
th  two  dependent  variables,  0  and  41  .  This  complication  is  normally 
noved  by  writing  KO^/az)  as  D(d0/dz)  ,  where  D  =  KCS'i'/a©)  is  defined 

the  soil  mositure  ditfusivity,  or  by  writing  90/'3t  as  C(9f/9t)  where 
=  90/9V  is  defined  as  the  soil  specific  moisture  capacity.  The  former  is 
ferred  to  in  the  following  as  the  Richards  equation  in  diffusivity  formula- 
on,  while  the  latter  is  referred  to  as  the  equation  in  pressure  head 
rmulation.  Either  formulation  may  be  used  for  most  applications,  but  if 
nded  water  on  the  surface  is  considered,  the  latter  form  is  required. 

8.  The  second  complication  is  that  the  coefficients  of  the  partial 
rivatives  in  either  formulation  are  strong  functions  of  the  dependent 
riables,  making  Equation  Cl  a  nonlinear  partial  differential  equation. 

is  complication  is  enriched  by  the  presence  of  hysteresis  and  by  nonlinear 
riation  of  the  partial  derivative  coefficients  over  a  range  of  5  to  7 
ders  of  magnitude  in  cases  of  interest.  The  nonlinearity  may  be  treated 

using  coefficients  at  each  time  step  which  are  derived  from  estimates  of 
ie  dependent  variables  at  the  next  time  step,  or  from  coefficients  appro¬ 
bate  to  the  dependent  variable  value  at  the  present  time,  or  from  some 
erage  value,  or  from  a  preliminary  calculation  of  the  new  dependent 
triable  value  which  is  used  to  calculate  coefficients  for  use  in  calculating 
corrected  dependent  variable  value,  or  by  iteration  in  which  new  coeffi- 
tnts  are  calculated  from  preliminary  calculations  of  the  dependent  variable 
ilues  until  two  iterations  result  in  values  which  are  sufficiently  close, 
ie  method  selected  by  each  modeler  is  generally  noted  in  the  following 
idel  discussions. 

9.  Where  the  model  is  formulated  in  such  a  way  that  hysteresis 
innot  be  ignored,  model  algorithms  generally  test  the  direction  of  soil 
listure  change  to  select  from  a  set  of  drying,  wetting  and  scanning  curves 
ie  one  which  is  appropriate  for  a  given  moisture  content  and  direction  of 
lange.  This  curve  is  then  used  in  a  nonlinearity  treatment  to  derive  the 
oper  coefficient  value  (or  value  of  matric  potential  from  moisture  content 
id  vice  versa)  . 
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10.  The  third  complication  arises  because,  as  a  partial  differential 
quation.  Equation  Cl  is  valid  only  for  infinitesimal  changes  in  time  and 
pace.  Digital  computers  are  incapable  of  infinitesimal  increments,  even  if 
ime  or  monetary  resources  are  available;  thus  selections  of  increments  of 
he  time  and  space  variables  for  realistically  accurate  computations  are 
equired.  In  most  of  the  models,  the  space  variable  increments  are  fixed, 
nd  the  time  increment  is  adjusted  by  the  program  to  meet  some  specifica- 
ion,  often  the  magnitude  of  change  of  a  dependent  variable  during  a  single 
ime  step.  In  some  models,  the  space  variable  increments  are  also  controlled 
hese  controls  enable  simulations  of  good  accuracy  without  excessive  expendi- 
ure  of  resources.  This  accuracy-related  control  of  increment  size  is  not 
lirectly  related  to  step  size  controls  for  stability  of  the  numerical  scheme 
rhich  are  discussed  below. 

11.  The  solution  of  Equation  Cl  by  numerical  methods  also  involves 
several  selections  among  alternate  formulations.  With  respect  to  time,  the 
ralue  of  the  dependent  variable  may  be  written  in  terms  of  values  known  from 
:he  previous  time  step,  or  in  terms  of  unknown  values  at  the  time  step  for 
rhich  the  calculation  is  being  made.  In  the  former  case,  the  unknown  value 
.s  an  explicit  function  of  known  values,  while  in  the  latter  case,  the 
mknown  value  is  an  implicit  function  of  itself  and  other  unknown  values, 
rhese  formulations  are  appropriately  identified  as  explicit  and  implicit, 
respectively.  Explicit  methods  are  the  most  straightforward  approach  and 
result  in  directly  solvable  expressions  for  each  time-space  grid  point,  but 
require  small  time  increments  to  avoid  instability.  On  the  other  hand. 
Implicit  methods  result  in  a  system  of  equations  which  must  be  solved 
simultaneously  for  all  grid  points  in  space  for  each  time,  although  they 

ire  absolutely  stable  in  linearized  systems.  In  addition  to  absolute 
stability,  the  implicit  methods  reward  the  greater  programming  effort 
required  by  actually  requiring  fewer  operations  to  resolve  the  entire  spatial 
»rid  each  time  step.  They  are  thus  doubly  efficient  in  comparison  t 
ixplicit  methods  by  permitting  larger  time  steps  (bounded  by  accuracy 
requirements  only)  and  simulating  each  step  more  rapidly. 

12.  Several  of  the  models  use  a  matrix  formulation  of  the  system  of 
simultaneous  equations  which  results  in  a  tridiagonal  coefficient  matrix 


the  implicit  formulation  in  one  dimension.  This  matrix  equation  may  be 
iently  solved  using  a  method  reported  by  Richtmyer  (1957).  Iteration 
ds  have  also  been  used  when  the  coefficients  were  not  linearized. 

13.  For  spatial  derivatives,  the  finite  difference  approximations  are 
ally  written  in  terms  of  the  values  of  the  dependent  variable  at 

ent  grid  points  on  either  side  of  a  point,  or  above  and  below  it. 

are  known  as  central  differences.  Central  differences  cannot  be  used 
e  boundaries  of  a  region  unless  imaginary  grid  points  are  introduced 
de  the  boundary;  thus  one-sided  spatial  differences  are  generally  used 
e  boundaries.  Boundary  conditions,  such  as  constant  matric  potential, 
cified  gradient  of  matric  potential,  a  no-flow  condition,  are  specified 
ie  model  applications  to  suit  the  physical  process  under  consideration, 
iriginal  papers  must  be  consulted  for  such  detail.  Further  general 
ission  is  provided  in  the  above-cited  references. 

14.  L.  F.  Richardson  has  been  credited  with  being  the  first  investi- 
■  to  seriously  apply  finite  difference  approximations  to  partial 
irential  equations  for  the  solution  of  flow  problems.  He  applied  the 

>d  to  a  problem  of  seepage  through  an  earth  dam  prior  to  1910.  He  is 
noted  for  his  use  of  numerical  methods  for  weather  forecasting  during 
following  decade.  However,  before  the  advent  of  electronic  digital 
iters,  such  pioneering  effort  was  severely  limited. 

15.  A.  Klute  (1952a  and  b)  is  generally  credited  with  developing  the 
:  numerical  model  for  solution  of  the  soil  moisture  flow  equation  by 
;al  computer.  In  keeping  with  a  suggestion  sketched  out  by  Childs  and 
Ls-George  (1950),  Klute  formulated  the  moisture  flow  equation  for 
netric  moisture  content  with  a  diffusivity  coefficient  as 


;  x  is  horizontal  distance  and  the  other  terms  are  defined  for  Equa- 
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Cl.  A  Boltzmann  transformation,  y  =  xt  ,  was  used  to  rewrite  the 
Lnear  partial  differential  equation  as  a  nonlinear  ordinary  equation  in 
Joltzmann  time-distance  variable,  y  .  An  iteration  scheme  was  then 
led  to  solve  for  the  soil  moisture  content  as  a  function  of  y  .  The 
L  simulation  of  flow  into  an  horizontal  soil  column  successfully 
iced  a  sharp  spatial  gradient  at  the  leading  edge  of  the  wetting  column 


:ing  front),  in  contrast  to  earlier  theoretical  models  of  the  phenomenon, 
mgh  the  front  appeared  to  move  too  rapidly,  and  weakened  with  time, 

»  defects  were  attributed  to  insufficient  input  data  leading  to  an 
iresis  problem.  Despite  these  defects,  a  significant  advance  had  been 


16.  Bruce  and  Klute  (1956)  undertook  the  measurement  of  diffusivity 

I  the  above  model.  They  were  able  to  calculate  moisture  diffusivity  as 
iction  of  soil  moisture  content  by  using  measured  values  of  moisture 
ent  against  distance  in  an  experimental  soil  column. 

17.  Ashcroft  et  al.  (1962)  published  a  model  for  the  solution  of 

e's  diffusion  equation  in  volumetric  moisture  content,  0  .  They  used  a 
y  implicit  technique.  The  equations  were  linearized  by  use  of  diffu- 
ty  values  for  each  time  step  calculated  by  the  method  of  Bruce  and  Klute 
6)  with  an  estimated  0  .  The  estimated  value  was  derived  by  adding  a 
1  increment  to  the  previous  0  value.  The  implicit  method  resulted  in 
stem  of  simultaneous  equations  which  were  solved  by  a  Gaussian  elimina- 
scheme.  The  model  was  applied  to  horizontal  diffusion  of  water  into  a 
geneous,  semi- inf inite  medium  with  good  results  in  comparison  to  experi- 
al  data,  and  to  analytical  solutions  utilizing  the  Boltzmann  transforma- 
.  While  this  model  has  not  had  a  significant  impact  on  succeeding  work, 
paper  is  notable  because  the  authors  included  rationale  and  discussion 
hoices  made  for  the  model  which  are  of  benefit  to  understanding  their 
and  the  work  of  many  others  who  have  been  quite  delinquent  in  this 
rd . 

18.  Whisler  and  Klute  (1965)  retained  the  iterative  characteristics 
he  earlier  Klute  model  in  a  new  formulation  of  the  problem,  but  their 
1  was  revised  in  many  important  respects.  The  revised  model  was  form¬ 
ed  to  calculate  the  time  rate  of  change  of  pressure  head,  h(=4')  , 

er  than  volumetric  moisture  content,  0  ,  as 


„/t_  ^  3h  9  ,  9h "1  .  9K(h,z) 

C(h-i!)  Si  ‘  S7  3iJ  +  -k-2 


(C3) 


e  all  symbols  have  been  defined  above.  This  equation  may  be  solved  for 
ral  vertical  flow  problems  with  appropriate  initial  and  boundary  condi- 
s,  and  was  used  in  several  of  the  following  applications  with  little 
f ication. 
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19.  The  Whisler-Klute  model  utilizes  nondimens ional  variables  and 
cients,  such  as  nondimensional  pressure  head,  given  by  $  =  h/L  , 

L  is  the  length  of  a  soil  column  under  study,  in  formulation  of  the 
ons  to  be  finite  differenced  and  solved.  Implicit  finite  differencing 
d  for  both  the  partial  derivative  approximations  and  the  coefficient 
imations  at  each  modeled  time  step.  Since  the  values  of  both  dimen- 
ss  pressure  head  and  the  equation  coefficients  are  thus  formulated  in 
of  the  values  at  adjacent  spatial  points,  a  system  of  nonlinear, 
aneous  equations  is  generated  for  the  values  at  the  given  time  step, 
et  of  algebraic  equations  cannot  be  solved  by  direct  elimination ,  or 
deterministic  techniques,  as  neither  the  coefficients  nor  the  pressure 
alues  are  known,  but  they  are  both  present  as  product  factors  in  the 
of  the  equation. 

20.  The  Whisler-Klute  model  system  of  algebraic,  nonlinear,  simulta- 
equations  is  solved  for  the  dimensionless  pressure  head  and  the 
.cients  by  iteration.  In  this  procedure,  values  for  the  coefficients 
'iven  time  are  first  approximated  by  calculation  of  their  values  using 
jad  values  from  the  prior  time  step.  These  approximate  coefficients 
len  used  as  constant  coefficients  in  the  set  of  simultaneous  equations, 
are  solved  for  pressure  head  at  each  spatial  grid  point  by  an 

:ified  technique.  The  improved  head  estimates  were  then  used  to  obtain 
iproximations  of  the  coefficients,  and  the  process  continued  until  the 
red  head  estimates  were  within  a  specified  amount  of  the  head  estimates 
:o  obtain  them.  The  final  set  of  coefficient  and  head  values  were 
led  after  this  convergence  and  used  to  calculate  the  coefficient  and 
values  for  the  succeeding  time  step  by  repeating  the  iteration 
5S.  The  entire  time  and  space  dependence  of  pressure  head  was  thus 
/ed . 

21.  The  theories  of  Childs  and  Collis-George  (1950)  and  Millington 
ilrk  (1961)  were  used  for  computer  calculation  of  the  necessary  varia- 
)f  hydraulic  conductivity  versus  head  from  curves  of  volumetric 

ire  content  versus  head.  The  latter  curves  were  arbitrarily  dravn  by 
ithors  for  typical  soil  types  as  relationships  which  appeared  reason- 
to  them.  The  combined  typical  soil  -dimensionless  analysis  approach  was 
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used  to  enhance  the  generality  of  the  results  of  model  simulations,  as 
specific  cases  of  soil  type  and  flow  geometry  with  matching  dimensionless 
variables  produce  identical  distributions  of  pressure  head  in  time  and 
space. 

22.  Whisler  and  Klute  (1965)  apply  the  model  to  the  analysis  of 
infiltration  into  a  column  of  soil  which  has  been  drained  to  equilibrium 
with  a  water  table  from  a  saturated  condition.  Whisler,  Klute,  and  Millington 
(1968)  use  a  similar  model  to  analyze  steady-state  evaporation  from  a  soil 
column,  including  water  uptake  and  redistribution  by  roots.  In  this 
application,  the  iteration  is  performed  for  the  single  time  representing  the 
steady-state  condition,  but  with  the  additional  complication  of  a  space 
varying  source  term  due  to  the  roots.  A  macroscopic  approach  to  root  water 
uptake  is  used,  where  only  bulk  characteristics  are  treated,  contrary  to  the 
approaches  of  Gardner  (1960)  and  Cowan  (1965) . 

23.  Whisler  and  Klute  (1969)  used  their  1965  model  to  study  addi¬ 
tional  cases  of  infiltration,  including  infiltration  from  rainfall  in 
addition  to  ponded  water  on  the  surface,  and  layered  soils.  Infiltration 
and  outflow  rates  were  also  calculated  with  the  19b9  model,  and  cumulative 
infiltration  was  included.  However,  the  finite  differencing  and  the 
iterative  model  solution  algorithm  were  essentially  unchanged. 

24.  Whisler  and  Watson  (1969)  modified  the  1965  Whisler-Klute  model 
to  include  provision  for  a  variable  time  step  in  a  simulation  of  infiltra¬ 
tion  into  a  draining  porous  medium.  The  variable  time  step  enabled  optimum 
computational  efficiency  while  ensuring  stable,  accurate  calculations. 

Their  time  step  was  not  set  arbitrarily,  however,  but  rather  was  adjusted  on 
the  basis  of  the  magnitude  of  change  in  pressure  head  for  the  previous  time 
step.  If  the  head  change  was  too  large,  the  step  was  repeated  with  a 
smaller  time  increment. 

25.  Watson  and  Whisler  (1972)  used  the  Whisler-Watson  1969  model  to 
treat  drainage  of  heterogeneous  media  from  saturation  with  a  water  table 
below,  but  no  surface  input  flux.  The  media  were  assumed  to  exhibit  scale 
heterogeneity  (see  Appendix  A).  While  field  heterogeneity  required  modi¬ 
fications  in  the  coefficient  expressions  of  the  mathematical  model,  the 
method  of  solution  of  the  equations  during  simulation  was  little  changed 
from  the  earlier  paper.  Whisler,  Watson,  and  Perrens  (1972)  used  the 


Whisler-Waston  1969  model  to  simulate  Infiltration  of  ponded  water  into  two 
soil  columns  with  different  heterogeneity  distributions  and  one  homogeneous 
column  for  comparison.  Hysteresis  was  avoided  in  this  study  by  specifying 
an  initial  moisture  content  that  was  constant  with  depth. 

26.  The  1965  Whisler-Klute  version  of  the  model  was  used  by  Klute  and 
Heermann  (1974)  to  simulate  soil-water  profile  response  to  periodic  surface 
boundary  conditions,  with  and  without  a  water  table.  A  high  degree  of 
harmonic  distortion  was  simulated  as  the  periodic  disturbance  propagated 
into  the  soil,  due  to  the  nonlinearity  of  the  soil-water  flow  system. 

27.  Bruce,  Thomas,  and  Whisler  (1976)  used  the  1969  Whisler-Watson 
model  version  to  predict  infiltration  into  layered  field  soils  to  study  the 
effects  of  various  distributions  of  hydraulic  characteristics  with  depth. 

The  soils  used  in  this  study  were  characteristic  of  field  soils,  contrary  to 
the  more  general  soil  types  studied  in  earlier  model  applications. 

28.  This  sequence  of  papers  indicates  the  range  of  problems  which  may 
be  treated  with  a  relatively  consistent  modeling  approach.  Similar  sequ¬ 
ences  of  work  by  several  authors  will  be  outlined  in  following  paragraphs 
for  several  different  approaches. 

29.  Philip  (1955,  1957a)  published  numerical  solution  models  for  the 
horizontal  diffusion  equation  and  for  vertical  flow  (summarized  in  Philip, 
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1957b).  His  solutions  also  utilized  the  Boltzmann  transformation,  <J>  =  xt  , 
to  reduce  the  partial  differential  equations  to  ordinary  differential  equa¬ 
tions,  as  did  Klute  (1952a),  but  Philip  then  made  volumetric  moisture 
content,  0  ,  the  independent  variable  and  formed  an  integro-differential 
equation  in  <J>  and  0  to  be  solved.  This  equation  was  solved  by  foward 
integration  with  an  initial  condition  derived  by  iteration.  The  vertical 
flow  case  required  an  approximation  of  distance  and  an  initial  estimate  by  a 
solution  to  the  equation  without  gravity.  A  series  of  partial  solutions  was 
formed  which  converged  to  the  final  solution.  As  this  approach  had  little 
Influence  on  succeeding  numerical  modeling  effort,  it  will  not  be  discussed 
further. 

30.  Another  early  soil  moisture  modeling  effort  which  did  not  exert 
great  influence  on  following  work,  but  is  of  historical  interest,  was 
published  by  Day  and  Luthin  (1956).  Their  numerical  solution  was  applied  to 
a  drainage  problem  for  which  they  had  obtained  laboratory  experiment  data. 

For  each  time  step  (treated  as  the  dependent  variable)  a  vertical  distribution 


of  pressure  head  was  estimated.  The  estimates  were  used  to  specify  K  and 
6  for  each  level,  and  to  calculate  the  water  loss  required  by  the  6 
changes  from  the  preceding  time  step.  An  explicit  centered  difference 
approximation  of  the  flow  equation  with  explicit  linearization  of  K  using 
previous  estimates  was  used  at  each  level  to  correct  the  head  value.  The 
cycle  was  iterated  until  the  distribution  of  head  values  did  not  change 
significantly,  and  the  time  lapse  consistent  with  the  top  head  value  changed 
(not  iterated)  was  calculated  from  the  total  desorption  in  the  column  and 
the  calculated  outflow  rate.  Steeping  through  time  was  accomplished  by  a 
new  (lower)  head  estimate  at  the  surface,  repeating  the  above  procedure. 

31.  While  the  introduction  of  the  high-speed  digital  computer 
rendered  this  procedure  obsolete,  several  of  the  more  sophisticated  models 
still  used  elements  of  the  method,  enhancing  its  historical  significance. 

The  model  predicted  higher  drainage  rates  than  measured,  and  therefore  a 
more  rapid  decrease  of  pressure  head,  although  it  reproduced  some  vertical 
gradients  caused  by  uneven  packing  of  the  soil  column  moderately  well. 
Departures  were  attributed  by  the  authors  to  both  uncertainties  in  the  K  -  0 
relationship,  determined  experimentally,  and  the  large  time  steps  dictated 
by  available  computational  means. 

32.  Hanks  and  Bowers  (1962)  were  the  first  to  combine  the  Crank- 
Nicholson  finite  difference  scheme,  a  rapid  algorithm  for  tridiagonal  matrix 
solution,  and  the  computing  power  of  a  digital  computer  in  a  numerical  model 
for  solution  of  the  soil  moisture  flow  equation.  Klute  (1952b)  had  suggested 
that  the  method  of  Crank  and  Nicholson  (1947)  might  be  adaptable  to  the  soil 
moisture  equation  for  vertical  flow,  but  after  working  out  the  necessary 
procedure,  Philip  (1957a)  found  considerable  labor  would  be  required  even 
for  limited  accuracy.  However,  by  programming  the  solution  for  processing 
by  an  IBM  650  computer  with  an  algorithm  for  solution  of  the  tridiagonal 
matrix  generated  in  the  Crank-Nicholson  method  that  was  published  by 
Richtmyer  (1957) ,  both  accuracy  and  labor  problems  were  brought  under 
control.  While  this  model  has  been  extended  in  many  ways  since  its  incep¬ 
tion,  and  while  other  numerical  methods  have  now  been  used  with  good 
results,  even  this  early  numerical  model  was  able  to  reproduce  moisture 
profiles  predicted  by  state-of-the-art  analytical  methods  with  good  accuracy; 
then  it  was  used  to  produce  profiles  for  layered  soils  which  were  well 
beyond  any  other  computational  capability. 


33.  The  Hanks-Bowers  numerical  model 
equation  in  one  dimension  in  the  form 

» i  i-  4 

at  c  az  y 

where 

h  ■■  pressure  or  tension  head 
t  *  time 

C  ■  specific  moisture  capacity  (C  * 
z  *  vertical  distance 
K  =  conductivity 
H  *  hydraulic  head  (H  =  h  +  z) 

0  *  volumetric  moisture  content 

34.  Specifics  of  the  finite  differencing  of  Equation  C4  are  published 
in  Hanks  and  Bowers  (1962),  which  should  be  consulted  for  details.  Basic¬ 
ally,  the  method  steps  forward  from  the  previous  to  the  present  time  using 
pressure  head  values  from  both  previous  and  present  time  levels.  It  is 
therefore  an  implicit  method  (explicit  methods  use  only  past  values  to  step 
to  the  present) .  A  set  of  simultaneous  equations  results  from  an  implicit 
method  at  each  time  step,  as  unknown  values  of  pressure  head  at  adjacent 
spatial  points  are  required  to  calculate  the  pressure  head  at  each  point. 
While  explicit  methods  are  conceptually  more  direct,  the  implicit  method  has 
better  stability  characteristics.  Further,  the  system  of  simultaneous 
equations  can  be  written  in  matrix  form  with  a  tridiagonal  coefficient 
matrix;  that  is,  only  the  main  diagonal  and  the  one  just  above  and  below  it 
have  nonzero  entries.  Using  efficient  algorithms,  the  system  of  simultaneous 
equations  can  actually  be  solved  with  fewer  operations  per  time  step  than 
are  required  for  the  apparently  simpler  explicit  method  (Remson,  Hornberger, 
and  Molz,  1971) . 

35.  The  need  for  iteration  to  resolve  both  coefficients  and  head 


solves  the  soil  moisture  flow 


(C4) 


ae/9h) 


values  was  eliminated  in  the  Hanks-Bowers  model  by  writing  C  and  K  in 
explicit  form  at  half-time  steps.  A  pseudo- imp lie it  form  is  actually  used 
for  C  ,  where  an  extrapolation  of  the  present  value  of  0  is  made  from 
prior  values,  and  C  is  calculated  using  this  estimated  value.  Superior 
results  were  obtained  when  K  was  written  as  a  somewhat  complicated  func¬ 
tion  of  soil  moisture  diffusivity,  rather  than  simpler  forms  also  tried. 


Cll 


Both  K  and  C  were  considered  constant  over  a  time  interval,  therefore 
linearizing  the  flow  equation  and  enabling  the  system  of  simultaneous  equa¬ 
tions  to  be  solved  without  iteration. 

36.  The  time  interval  was  made  variable  in  order  to  achieve  both  good 
temporal  resolution  during  the  rapidly  changing  infiltration  periods  and 
computational  efficiency  during  later  periods  of  slow  change.  The  interval 
was  defined  as  the  time  for  a  specified  quantity  of  water  to  enter  the  soil, 
and  was  adjusted  after  each  time  step.  Cumulative  infiltration  was  computed 
from  the  change  of  9  ,  while  infiltration  rate  was  calculated  from  the 
pressure  head  gradient  at  the  surface  and  the  conductivity  in  the  first  soil 
layer . 

37.  The  depth  increments  were  defined  for  the  layered  soil  cases  such 
that  the  change  of  soil  type  occurred  at  a  half-space  level.  Presure  head 
was  required  to  be  continuous  across  the  boundary,  as  was  flow;  thus  volu¬ 
metric  moisture  content  changed  abruptly  at  the  boundary  of  the  differing 
soils.  The  conductivity  at  the  half-step  was  defined  as  an  average  of 
conductivities  above  and  below  the  boundary. 

38.  The  model  requires  initial  and  boundary  conditions  and  known 
relations  between  moisture  content  and  both  conductivity  and  diffusivity.  A 
series  of  computations  are  made  at  each  time  step,  and  the  model  is  stepped 
through  time  for  a  specified  period  to  effect  the  solution  desired. 

39.  The  model  was  used  by  Hanks  and  Bowers  (1962)  to  calculate 
infiltration  into  homogeneous  and  layered  media.  The  homogeneous  media 
results  for  two  soils  were  compared  to  solutions  for  the  same  soils  derived 
analytically  after  Scott  et  al.  (1962)  and  Philip  (1955).  In  these  compari¬ 
sons  the  same  restrictive  assumptions  required  by  the  analytical  methods 
were  applied  to  soil  characteristics  input  to  the  model.  Very  good  corre¬ 
spondence  between  analytical  and  numerical  results  was  found,  except  for 
somewhat  erratic  simulated  infiltration  rates  at  midtime  of  the  simulations. 
The  model  was  then  applied  to  layered  soil  cases  (fine  over  coarse  soil  and 
coarse  over  fine  soil).  The  results  were  reasonable,  but  data  for  compari¬ 
son  were  not  available. 

40.  Hanks  and  Bowers  (1963)  applied  the  model  to  evaluate  the 
influence  of  variations  in  the  dif fusivity-watet  content  relations  on 
infiltration,  while  Hanks  and  Gardner  (1965)  used  it  to  evaluate  the 
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influence  of  the  variations  on  surface  evaporation.  In  the  latter  case,  the 
model  was  modified  to  handle  surface  drying  due  to  evaporation.  Simulated 
model  results  were  found  to  be  consistent  with  other  research  results,  and 
indicated  that  accuracy  in  the  relation  between  diffusivity  and  water 
content  is  most  important  at  high  water  contents. 

41.  Jensen  and  Hanks  (1967)  applied  the  model  to  the  problem  of 
nonsteady-state  drainage  from  porous  media,  with  comparison  to  laboratory 
data.  The  method  of  determining  K  was  modified  for  this  application,  but 
the  model  was  otherwise  unchanged.  Very  good  correspondence  with  three 
laboratory  experiments  on  different  soils  was  achieved,  except  at  very  short 
times.  Some  discrepancy  between  data  and  simulation  could  also  be  attributed 
to  imperfections  in  the  soil  packing  for  the  experiment,  etc.;  however,  as 
these  details  were  unknown,  they  could  not  be  included  in  the  model  input. 

42.  The  complication  ot  hysteresis  in  the  pressure  head-water  content 
relationship  was  avoided  in  the  above  model  applications  by  starting  from  a 
moisture  content  constant  with  depth  and  only  wetting  (infiltration)  or 
drying  (evaporation  and  drainage)  the  soil.  This  restriction  on  the  model 
was  removed  in  1969  (Hanks,  Klute,  and  Bresler,  1969),  and  the  model  was 
applied  to  a  complex  case  of  infiltration,  redistribution,  drainage,  and 
evaporation  from  the  soil.  While  the  K-8  relationship  was  assumed  to  be 
without  hysteresis  (generally  valid) ,  both  C-9  and  h-0  relations  were 
treated  as  hysteretic.  The  method  used  was  similar  to  that  of  Rubin  (1967). 
This  required  a  model  modification  to  evaluate  whether  wetting  or  drying  was 
occuring  during  a  time  step  and  whether  or  not  this  change  was  the  same  as 
the  prior  time  step.  The  appropriate  scanning  curve  of  the  hysteretic 
relationship  was  then  selected  to  derive  8  and  C  from  computed  values  of 
h  .  With  this  modification  and  specification  of  appropriate  initial  and 
boundary  conditions,  the  model  was  applied  to  the  complex  case.  Model 
simulation  results  were  compared  to  experimental  data  from  laboratory  soil 
columns  wetted  at  three  different  rates,  then  allowed  to  drain  with  and 
without  evaporation  from  the  surface  (Bresler,  Kemper,  and  Hanks,  1969). 
Evaporation  was  underestimated  somewhat,  while  drainage  during  evaporation 
was  somewhat  overestimated,  but  the  results  of  the  comparison  were  generally 
very  good. 


C13 


43.  The  1969  version  of  the  model  was  used  by  Bresler  and  Hanks 
(1969),  with  an  added  algorithm  for  treatment  of  salt,  to  estimate  simul¬ 
taneous  movement  of  water  and  salt  in  unsaturated  soils.  The  same  problem 
was  treated  by  Warrick,  Biggar,  and  Nielsen  (1971),  using  the  1962  Hanks- 
Bowers  model  with  a  different  salt  algorithm.  In  each  case,  field  data  or 
laboratory  data  was  used  for  comparison  with  simulations  by  the  model.  (In 
a  subsequent  paper,  Reichardt,  Neilsen,  and  Biggar  (1972b)  used  a  new 
explicit  model,  with  C  and  K  definitions  similar  to  those  of  Hanks  and 
Bowers,  to  study  horizontal  infiltration  into  layered  soils.)  The  1969 
version  of  the  model  was  also  used  by  Bresler  (1973),  although  with  an 
extensively  revised  algorithm  for  treatment  of  salt  movement  without 
numerical  dispersion.  Again,  these  modeling  efforts  were  compared  to 
laboratory  and  field  data  with  good  correlation. 

44.  While  modeling  of  soil  moisture  under  bare  fields  is  of  value  in 
many  applications  and  aids  the  understanding  of  soil  water  movement  under 
those  conditions,  it  is  not  sufficient  for  agriculture  or  for  more  general 
concerns  of  traff icability  and  hydrology.  As  noted  in  the  discussion  of 
evapotranspiration  (Appendix  B) ,  soil  water  depletion  via  evaporation  and 
transpiration  is  a  complex  interrelated  set  of  processes,  which  include 
plant  influences  in  both  the  soil  and  the  atmospheric  boundary  layer. 
Effective  soil  moisture  modeling  must  include  the  influences  of  plants, 
either  implicitly  in  the  determination  of  parameters,  or  explicitly. 

45.  Nimah  and  Hanks  (1973a  and  b)  applied  the  Hanks-Bowers  1962  model 
to  the  problem  of  soil  moisture  prediction  under  an  alfalfa  crop.  The  model 
was  expanded  for  this  application  by  including  a  root  extraction  function 
for  soil  water,  partitioning  potential  evapotranspiration  into  potential 
evaporation  and  potential  transpiration  to  treat  the  differing  depletion 
mechanisms  and  surface  boundary  conditions,  and  including  an  iteration 
procedure  to  solve  for  the  changing  surface  boundary  conditions.  Exchange 
of  moisture  with  a  water  table  was  included  in  the  model,  but  hysteresis 
effects  were  ignored.  Potential  evapotranspiration  was  determined  with 


a  modified  Penman  approach.  The  numerical  procedure  was  also  modified 
to  include  a  variable  depth  increment  to  improve  efficient  resolution  of 
gradients,  while  the  variable  time  increment  was  based  on  the  total  change 
of  water  content  in  the  column  for  this  model  application.  Brief  compari¬ 
son  with  field  data  is  made  in  the  first  of  the  two  papers,  while  more 
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extensive  comparison  is  made  in  Nimah  and  Hanks  (1973b)  against  data  from  a 
field  experiment.  Model  simulations  compared  well  in  both  time  and  depth, 
except  for  immediately  after  irrigation  or  rainfall. 

46.  Feddes,  Bresler,  and  Neuman  (1974)  used  a  modified  form  of  tha 
1973  Nimah-Hanks  model  to  simulate  water  balances  in  comparison  to  field 
data  on  root  uptake  of  soil  water  by  red  cabbage.  Their  modifications 
included  the  surface  layer  under  dry  surface  conditions,  rather  than  a  fixed 
air-dry  water  content,  use  of  a  potential  evapotranspiration  model  which 
included  crop  structure  characteristics,  and  use  of  a  revised  root  extrac¬ 
tion  function.  They  found  that  the  changed  root  extraction  function  and  the 
use  of  an  equilibrium  surface  moisture  content  made  little  difference  on 
simulated  values  of  evaporation,  transpiration,  and  soil  moisture,  but  that 
the  changes  of  the  potential  evapotranspiration  estimation  were  of 
significant  value. 

47.  Childs  and  Hanks  (1975)  added  diurnal  variation  of  potential 
evapotranspiration  and  seasonally  varying  partitioning  of  potential  evapora¬ 
tion  and  potential  transpiration  on  the  basis  of  crop  development  to  the 
Nimah-Hanks  model.  The  approach  of  the  latter  modification  was  extended  by 
Childs,  Gilley,  and  Splinter  (1977)  in  a  simple  model  which  included  the 
above-ground  plant  functions,  as  well  as  root  extraction.  The  Childs-Hanks 
model  was  used  by  Watts  and  Hanks  (1978)  with  the  Bresler  (1973)  formulation 
of  solute  transfer  to  study  soil-water-nitrogen  interactions.  Comparison 
with  field  data  was  made  in  each  case  with  good  correlation,  but  neither 
perfect  nor  completely  consistent  correlation,  as  in  the  other  cases. 
Tillotson  et  al.  (1980)  report  the  current  state  of  the  model  of  soil 
moisture,  plant  root  interaction,  and  solute  flow.  This  publication 
provides  a  more  complete  exposition  of  the  model  than  do  the  journal  papers. 

48.  Liakopoulos  (1966)  used  a  revised  formulation  of  Equation  C4  to 
model  unsteady,  unsaturated  soil  moisture  flow.  He  expanded  hydraulic 
head  H  *  h  +  z  ,  wrote  out  the  derivative  on  the  right  side  of  Equation  C4, 
used  the  standard  definition  of  soil  moisture  diffusivity,  D  =  K/C  = 
K(9h/30)  ,  defined  B  =  3K/99  ,  and  derived  the  equation 


where 


h  =  pressure  head 
t  =  time 

z  =  vertical  coordinate 
3  =  volumetric  moisture  content 

The  time  derivative  was  written  as  a  finite  forward  difference,  while  the 
second  derivative  was  written  as  an  implicit  centered  difference,  and  the 
first  derivatives  of  the  second  term  on  the  right  were  written  as  explicit 
centered  differences.  The  equation  as  nonlinearity  was  removed  by  writing 
the  coefficients  D  and  B  in  explicit  form,  evaluated  at  the  beginning  of 
the  time  step. 

49.  The  implicit  finite  difference  formulation  results  in  a  system  of 
linear  equations  to  be  solved  for  the  values  of  pressure  head  at  all  grid 
locations  each  time  step.  This  system  was  solved  by  Gaussian  elimination 
using  the  same  matrix  formulation  as  used  by  Hanks  and  Bowers  (1962)  and 
others,  but  Liakopoulos  fully  described  the  solution  scheme,  which  has 
otherwise  been  covered  only  in  textbooks.  The  derivation  of  the  tridiagonal 
coefficient  matrix  and  the  recursive  method  for  determination  of  coeffi¬ 
cients  of  a  reduced  matrix  are  clearly  presented,  and  are  recommended 
reading,  as  this  scheme  is  quite  frequently  used  in  soil  moisture  modeling. 

50.  Liakopoulos  claims  satisfactory  comparison  of  results  from  the 
model  with  actual  experimental  data,  but  presents  only  theoretical  results 
for  gravity  drainage,  evaporation,  infiltration,  and  capillary  rise. 

Initial  conditions  were  formulated  in  each  case  to  avoid  hysteresis. 

51.  Staple  (1966)  solved  the  moisture  flow  equation  using  explicit 
finite  differences  in  both  the  original  form  derived  by  Richards  (1931), 
namely , 


where  V  =  pressure  head  or  matric  potential  (h  above)  and  the  other  symbols 
are  defined  as  above;  and  in  £he  form  derived  by  Klute  (1952a),  namely. 


where  D  =  K(94'/30)  *  K/C  ,  as  above.  Staple  used  Equation  C6  for  the  upper 
profile  ot  the  soil  during  drying,  and  he  used  Equation  C7  during  infiltra¬ 
tion  of  the  soil  column  and  wetting  of  the  lower  soil  column.  Staple  (1966) 
cites  an  early  work.  Staple  and  Lehane  (1954),  in  which  Equation  C6  was 
solved  by  explicit  finite  difference.  Historically,  this  would  place  his 
work  just  after  that  of  Klute.  The  1954  paper  is  seldom  cited  by  others, 
however,  and  authors  have  generally  chosen  to  eliminate  either  0  or  'f(h) 
from  the  Richards  equation  before  proceeding  to  finite  different ing  and 
numerical  solution. 

52.  Staple  (1966)  used  average  values  of  conductivity,  K  ,  over 

adjacent  grid  points  and  the  method  of  Staple  and  Lehane  (1954)  and  Hanks 
and  Bowers  (1962)  to  derive  mean  values  of  diffusivity,  D  ,  for  use  in  the 
solution.  Laboratory  data  was  used  for  the  K-0  ,  K-^  ,  and  D-0  relations 

required,  obtained  from  earlier  work  with  emphasis  on  hysteresis  (Staple, 

1962  and  1965) .  Tabulated  input  data  were  used  in  the  computer  program. 

53.  Staple  (1966)  applied  his  model  to  simulation  of  moisture 
profiles  during  infiltration  and  subsequent  redistribution  in  a  soil  column. 
The  results  were  found  to  be  of  the  right  magnitude  in  comparison  to  related 
measurements.  Staple  (1969)  applied  the  same  model  to  measured  laboratory 
data  (with  an  improved  algorithm  for  treatment  of  the  transition  from 
wetting  to  drying)  and  found  the  agreement  satisfactory.  In  this  latter 
application.  Staple  used  implicit  formulation  of  the  flow  equation  in  some 
of  the  tests,  since  it  was  faster. 

54.  Rubin  (1969)  (publication  of  1966  Symposium  paper  in  Rijtema 
and  Wassink,  1969)  presented  another  approach  to  numerical  solution  of  the 
moisture  flow  equation.  He  reduced  the  Richards  equation  in  two  variables, 

0  and  f  ,  to  an  equation  in  a  single  dependent  variable,  v  ,  via  a 
Kirchhoff  transformation: 


v  -  v(H) 


K(h)dh 


max 


(C8) 


where 


and  where  both  H  and  h  represent  the  pressure  head  here,  with 
and  the  upper  and  lower  bounds  of  pressure  heads  in  the  medium  under 

consideration.  Assuming  the  inverse  of  v(H)  exists  and  3H/3v  =  V/K(H)  , 
the  Richards  equation  may  be  written  in  the  single  dependent  variable,  v  ,  as 

2 

Y(v)  it  =  Yi  ~  z(v)  37  (cl0) 

Z 

where 

Y(v)  =  C/K(H(v) ) 

Z(v)  =  3 In  K(H)/3H 
C  =  30/3H 

(See  Rubin's  paper  or  a  general  text  discussion  of  the  Kirchhoff  trans¬ 
formation  for  more  detail.) 

55.  Rubin's  justification  of  use  of  the  Kirchhoff  transformation 
instead  of  the  pressure  head  formulation  is  not  confirmed  by  the  comparison 
evalations  performed  by  Havercamp  et  al.  (1977),  but  his  justification  of 
its  use  instead  of  the  moisure  content  formulation  has  been  noted  elsewhere 
in  this  report.  Actually,  it  is  a  viable  alternate  approach  which  may  be 
preferred  under  certain  circumstances;  it  needs  no  further  justification. 

56.  The  numerical  solution  is  effected  by  means  of  the  Crank- 
Nicholson  implicit  scheme,  which  is  also  used  by  several  other  investi¬ 
gators.  Actually,  the  Crank-Nicholson  scheme  involves  an  averaging  of  both 
implicit  and  explicit  terms,  but  is  generally  referred  to  as  implicit 
without  qualification.  Nonlinearity  due  to  the  coefficients  Y  and  Z 

is  removed  by  an  explicit  determination  of  their  values  at  the  midpoint  of 
the  time  step  from  known  values  of  v  ,  Y  ,  and  Z  at  the  beginning  of  the 
time  step.  First,  v  is  calculated  for  the  half-step,  then  Y(v)  and 
Z(v)  are  computed.  K  is  also  required  in  the  finite  difference  expression 
for  3v/3z  ,  and  is  calculated  with  an  explicit  method  for  the  entire  time 
step.  The  specific  method  of  solving  the  resulting  system  of  linear 
simultaneous  equations  is  not  noted  in  the  paper,  but  is  likely  to  have 
been  the  efficient  tridiagonal  matrix  algorithm  in  Richtmyer  (1957),  used  by 
others. 

57.  The  program  uses  tabulated  values  of  Y(v)  and  Z(v)  derived 
from  functions  fitted  to  data  for  the  soils  under  consideration.  These 
functions  are  also  used  to  derive  h  output  values  from  the  computed 

v  values. 


58.  This  numerical  model  was  applied  to  calculation  of  the  soil 
joisture  profiles  during  preponding  infiltration  of  rainfall  and  infiltra- 
:ion  from  rain  ponds  on  the  surface.  The  results  were  found  both  reasonable 
ind  enlightening,  although  only  profiles  under  flooding  were  available  for 
:omparison. 

59.  Rubin  (1967)  used  a  different  model  for  analysis  of  postinf iltra- 
:ion  redistribution  of  soil  moisture.  In  this  model  Equation  C7  (in  6) 

/as  solved  using  implicit  finite  difference  formulation,  based  on  earlier 
/ork  by  Rubin  and  Steinhardt  (1963) .  The  equation  was  linearized  by  extra¬ 
polation  of  0  from  prior  time  steps  and  calculation  of  D  and  K  from 
:he  extrapolated  values.  The  particular  method  of  solving  the  system  of 
Linear  simultaneous  equations  is  not  noted. 

60.  This  paper  is  notable  for  the  method  of  hysteresis  treatment.  As 
the  K-0  and  D-0  relationships  are  functions  of  whether  drying  or  wetting 
is  occurring  and  the  moisture  content  at  the  time,  the  condition  and  direc¬ 
tion  of  change  of  soil  moisture  must  be  considered  in  calculation  of  K 

and  D  from  0  .  To  this  end,  a  second  grid  is  defined  with  a  one-to-one 
correspondence  with  the  soil  moisture  grid  in  which  this  information  is 
held.  The  decision  as  to  use  of  a  wetting,  drying,  or  scanning  curve  is 
made  on  the  basis  of  this  information  each  time  D  or  K  is  determined, 
and  when  moisture  contents  are  converted  to  pressure  head  values.  The 
required  information  is  simplified  by  the  assumption  that  once  drying  begins 
at  a  grid  point  under  redistribution  it  will  continue.  Thus  the  coded 
information  is:  use  wetting  relationship  or  use  dryi  ig  or  scanning  curve 
with  parameter  stored  in  second  grid  point.  Details  may  be  found  in  Rubin 
(1967). 

61.  Rubin  (1963)  formulated  a  numerical  model  for  use  in  two- 
dimensional  cases  in  unsaturated  and  partly  saturated  soils.  The  Kirchhoff 
transformation  used  in  Rubin  (1967)  was  used  to  model  a  case  of  two- 
dimensional  horizontal  infiltration,  while 


(Cll) 


was  used  for  a  falling-water-table,  ditch-drainage  case,  where  H  =  h  +  z  , 
C  =  30/f)h  ,  the  specific  water  capacity,  K  is  the  hydraulic  conductivity, 
and  t  ,  x  ,  and  z  are  independent  variables  of  time,  horizontal  dis- 
:ance,  and  vertical  distance,  respectively. 


Lth  set  (2)  data  supplemented  by  field  examination  of  a  single  boring 
mdomly  placed  in  the  nodal  area,  including  rooting  depth;  and  (4)  simula- 
Lons  were  performed  with  additional  data  available  through  field  and 
aboratory  measurements  of  the  soil  from  the  single  boring  of  set  (3). 

103.  They  found  simulation  of  the  water  table  depths  to  be  poor  for 
ata  sets  (l)-(3),  but  to  be  greatly  improved  using  data  set  (4).  This 
mprovement  proved  to  be  due,  primarily,  to  better  characterization  of  the 
ubsoil  properties,  however,  as  results  with  data  sets  (l)-(3)  were  greatly 
mproved  by  the  single  addition  of  subsoil  data  from  set  (4).  These  results 
ay  in  part  be  due  to  the  two-layer  model  used,  also,  as  was  noted  for  the 
esults  of  simulation  of  evapotranspiration.  In  this  latter  case,  the  model 
as  quite  sensitive  to  rooting  depth.  To  further  check  on  the  representa- 
iveness  of  the  results,  and  in  particular  to  check  on  the  data  set  (2) 

nd  data  set  (4)  differences  in  simulated  evapotranspiration,  several 
dditional  borings  were  made  in  the  five  nodal  areas.  It  was  found  that 
he  single  boring  used  for  data  set  (4)  was  generally  less  representative  of 
he  field  conditions  of  the  nodal  area  than  the  estimates  made  from  the  soil 
iurvey  information  for  data  set  (2).  It  was  thus  recommended  not  to  perform 
my  soil  borings  in  the  field  unless  several  were  possible. 

104.  It  is  notable  that  the  principal  difficulty  encountered  in  this 
simulation  effort  was  caused  by  subsoil  characteristics  which  were  poorly 
•epresented  by  data  available  from  the  soil  surveys,  and  even  more  poorly 
•epresented  by  a  single  soil  boring.  Present  computer  hardware  and  software 
sapabilities  for  simulation  of  soil  moisture  regimes  have  so  far  outstripped 
she  available  data,  that  even  this  greatly  simplified  model  is  limited  by 
ivailable  inputs. 


Mass  Balance  Modeling  of  Soil  Moisture 


105.  Mass  balance  models  of  soil  moisture  have  an  appeal  which  derives 
rom  their  simplicity  of  concept  and  generally  uncomplicated  application, 
n  such  a  model,  one  simply  accounts  for  the  accretion  and  depletion  from  a 
ioil  column  or  soil  layer  by  addition  or  subtration  of  an  appropriate 
[uantity  of  water,  and  the  moisture  volume  remaining  after  each  step  s  the 
iroper  moisture  content.  This  simplicity  is  very  misleading,  however,  as 


with  field  data  due  to  Jackson  (1972)  was  obtained,  however,  indicating  to 
the  authors  that  the  model  equations  are  reasonable  representations  of  the 
physical  processes  modeled.  One  obvious  problem  may  be  noted,  however. 

They  used  the  simplified  expressions  of  Clapp  and  Hornberger  (1978)  for  the 
relations  of  pressure  head  and  hydraulic  conductivity  to  volumetric  moisture 
content,  and  their  Figure  3-7  shows  significnat  departures  of  these  expres¬ 
sions  from  the  field  data  at  high  water  contents.  Since  Hanks  and  Bowers 
(1963),  Hanks  and  Gardner  (1965)  and  Van  Keulen  and  Hillel  (1974)  had 
clearly  demonstrated  that  accuracy  in  these  quantities  is  particularly 
important  at  high  water  contents,  some  of  the  model  error  is  likely  to  be 
from  this  source. 

101.  Finally,  a  numerical  model  for  soil  moisture  simulation  by  De  Laat 
(1976)  was  used  by  Bouma  et  al.  (1980a)  to  treat  the  effects  of  lowering 
water  tables  on  grass  production  using  soil  survey  identification  of  soils, 
and  by  Bouma  et  al.  (1980b)  to  simulate  regional  soil  moisture  regimes  using 
soil  survey  data.  The  model  of  De  Laat  (1976)  is  characterized  more  by  its 
design  intent  of  great  computational  efficiency  than  by  its  sophistication, 
being  a  pseudo-steady  state,  two-layer  model  using  Penman  method  estimation 
of  evapotranspiration,  but  this  application  is  of  great  interest  for 
predicting  traf f icability  in  inaccessible  areas. 

102.  The  former  paper  by  Bouma  et  al.  (1980a)  includes  mapping  via 
computer  of  simulated  soil  sensitivity  to  water  table  drawdown.  They 
emphasize  the  advantages  of  computer  storage  of  soil  survey  and  simulation 
information  which  may  be  mapped  by  the  computer  in  response  to  specific 
concerns,  in  comparison  to  the  traditional  mapping  of  interpretations.  The 
latter  paper  by  Bouma  et  al.  (1980b)  discusses  a  sequence  of  simulations  for 
soils  identified  by  a  soil  survey  map.  An  area  of  6  km  by  6  km  (3.7  miles 
by  3.7  miles)  was  gridded  for  computer  modeling,  and  data  were  presented  for 
five  selected  nodal  areas  of  25  hectares  (62  acres)  each.  Simulations 

were  performed  for  comparison  with  field  data  using  four  data  sets  of 
increasing  specificity,  namely:  (1)  Simulations  were  performed  for  each  of 
the  five  nodes  using  soil  physical  characteristics  reported  in  the  litera¬ 
ture  for  the  soil  type  most  common  in  the  6  km  by  6  km  area;  (2)  simulations 
were  performed  for  the  most  common  soil  type  in  each  25  ha  nodal  area  using 
data  from  sources  as  in  set  1;  (3)  simulations  were  performed  for  each  area 


simulated  transpiration  from  a  corn  plant  quite  well,  once  the  initially 
assumed  plant  resistance  was  adjusted.  The  difficulty  of  plant  resistance 
determination  was  noted  as  a  model  limitation.  This  model  is  not  unique  in 
this  regard. 

98.  Two  papers  treating  ground  water  flow  in  relation  to  flow  in  the 
unsaturated  zone  were  also  noted,  despite  an  extremely  thin  connection  with 
agricultural  science,  and  are  briefly  included  here.  There  is  a  large  body 
of  literature  related  to  watershed  modeling  which  has  not  been  included 
because  of  the  focus  of  this  work  on  agricultural  developments,  and  these 
two  papers  provide  a  "window"  into  the  region. 

99.  Green  et  al.  (1970)  modeled  the  movement  of  water  under  a  shallow 
pond  using  an  implicit-iterative  technique.  They  considered  movement  of 
both  air  and  water  in  a  two-phase  flow.  They  obtained  good  correlation 
between  experiments1  data  and  model  simulation,  but  only  after  some  modifica¬ 
tion  of  the  originally  assumed  values  of  the  porous  media  properties.  As  in 
so  many  other  cases,  the  numerical  solution  to  the  partial  differential 
equations  used  to  model  the  complex  water  movement  process  is  far  less 
uncertain  than  knowledge  of  the  actual  physical  conditions  in  a  field 
experiment.  Freeze  (1971)  extended  his  one-dimensonal  modleing  effort 
(Freeze,  1969)  to  three  dimensions.  The  resulting  equations  are  solved  by  a 
line  successive  overrelaxation  technique,  and  the  model  is  quite  general  in 
treatment  of  a  small  watershed. 

100.  A  program  considering  the  simultaneous  flow  of  heat  and  moisture 
in  soils,  reminiscent  of  the  work  of  Van  Bavel  and  Hillel  (1976),  was 
developed  by  Camillo  and  Schmugge  (1981)  for  use  in  conjunction  with  remote¬ 
sensing  techniques  for  soil  moisture.  They  formulated  the  equations  of 
heat,  liquid  water,  and  water  vapor  flow  in  terms  of  respective  diffusi- 
vities,  including  the  movement  of  water  vapor  through  the  atmospheric 
boundary  layer.  The  general  solution  was  accomplished  via  an  Adams-Bashford 
finite  difference  approach,  while  the  nonlinear  dependence  on  surface 
temperature,  created  by  terms  of  the  surface  energy  budget  formulation,  was 
treated  via  iteration.  This  publication  is  notable  for  the  completeness  of 
their  description  of  the  numerical  solution.  Simulated  results  compare  well 
to  analytic  and  quasi-analytic  solutions,  where  the  precise  boundary  condi¬ 
tions  and  soil  properties  are  identical  for  the  two  approaches,  but  problems 
again  develop  in  application  to  field  conditions.  Qualitative  correspondence 


The  Millington  and  Quirk  (1961)  method  was  used  to  calculate  conductivity 
and  diffusivity  from  field-measured  relations  between  moisture  content  and 
matric  potential.  Model  results  compared  favorably  with  data. 

95.  The  successive  overrelaxation  (SOR)  technique  was  used  by  Amerman 
(1976)  and  Reisenauer  (1963)  to  treat  two-dimensional  and  multidimensional 
soil  water  movement,  respectively,  while  the  finite  element  approach  was 
used  by  Neuman  and  Witherspoon  (1970);  Guymon,  Scott,  and  Herrmann  (1970); 
Neuman,  Feddes,  and  Bresler  (1975),  theory;  Feddes,  Neuman,  and  Bresler 
(1975),  field  application;  and  Parkes  and  O' Callaghan  (1980).  A  numerical 
method  using  the  flow  velocity  equation  was  used  with  success  by  Wind  and 
van  Doome  (1975)  ,  and  Richter  (1980) .  Each  of  these  methods  has  value  in 
special  applications,  such  as  use  of  the  finite  element  technique  for 
problems  with  complex  boundaries.  However,  they  have  not  been  in  the 
mainstream  of  development  and  will  not  be  discussed  further  here.  Remsom, 
Homberger,  and  Molz  (1971)  provide  a  thorough  discussion  of  the  SOR  and 
finite  element  methods,  which  should  be  consulted  for  details. 

96.  Two  additional  papers  relating  to  the  modeling  of  water  uptake  by 
roots  should  be  noted.  Feddes  et  al.  (1976)  use  an  implicit  finite  differ¬ 
ence  model  of  the  flow  equation  in  diffusivity  formulation  with  volumetric 
moisture  content  as  the  dependent  variable  and  an  added  root  uptake  source 
term.  The  method  of  linearization  of  the  equations  during  calculation  is 
not  specified,  but  appears  to  be  linear  averaging  of  the  coefficients.  The 
root  effectiveness  function  used  was  a  function  of  soil  moisture  content  and 
critical  moisture  contents  for  root  activity.  Root  growth  was  considered. 
While  cumulative  evaporation  and  transpiration  were  well  simulated  in  compari¬ 
son  to  field  data,  the  vertical  distribution  of  soil  water  content  was  not. 
This  simpler  model  did  compare  favorably  in  simulation  of  field  measurements 
with  the  more  complex  model  of  Feddes,  Bresler,  and  Neuman  (1974),  although 
the  two  models  frequently  produced  opposite  departures  from  the  data. 

97.  The  second  paper  was  published  by  Slack,  Haan,  and  Wells  (1977). 
They  used  a  microscopic  approach  to  evaluate  the  root  extraction  function  of 
depth,  then  a  macroscopic  model  to  solve  the  flow  equation  with  the  added 
sink  term.  An  implicit  method  was  used,  but  details  were  given  by  reference 
to  a  Ph.D.  dissertation  which  was  unavailable  for  this  review.  The  model 


values  by  Van  Bavel  and  Ahmed  (1976)  revealed  higher- than-normal  evapotrans- 
piration.  However,  as  they  had  modeled  a  sequence  of  fair-weather  June 
days,  the  excess  may  have  been  due  to  that  choice. 

92.  A  different  continuous  simulation  language,  DYNAMO  II  (Pugh, 

1970),  was  utilized  by  Hansen  (1975)  to  formulate  a  numerical  model  of  the 
water  state  and  transportation  in  the  soil-plant-atmosphere  system.  While 
an  appendix  contains  a  complete  model  listing,  no  discussion  of  the  relative 
merits  of  DYNAMO  II,  CSMP,  and  modeler-written  computer  algorithms  is 
included.  It  is  noted  that  the  basic  integration  follows  Euler's  methods, 
thus  explicit  methods  are  used.  The  model  addresses  the  problem  of  the 
state  and  flow  of  water  in  the  soil  and  growing  crop,  and  is  thus  time 
dependent.  It  was  partially  based  on  experimental  data,  but  discussion  of 
comparison  to  data  was  deferred  to  a  subsequent  paper  (not  reviewed) . 

93.  Several  other  modeling  approaches  to  the  problem  of  soil  moisture 
have  been  used,  but  do  not  fall  naturally  into  the  above  groups.  One  such 
approach  was  taken  by  Wang  and  Lakshminarayana  (1968)  to  simulate  water 
movement  in  nonhomogenous  soils.  The  Richards  equation  (Equation  Cl)  was 
used,  but  the  vertical  derivatives  were  first  written  after  the  fashion  of 
Liakopoulos  (Equation  C5)  with  further  consideration  of  the  spatial  varia¬ 
tion  of  hydraulic  conductivity,  due  to  the  application  to  heterogeneous 
soils.  Both  explicit  and  implicit  finite  difference  formulations  were  used, 
while  the  nonlinearity  of  the  equation  was  addressed  by  use  of  iteration. 
While  hysteresis  was  neglected  in  this  model  formulation,  calculated  and 
field -measured  data  compare  favorably  for  vertical  drainage  and  infiltration. 

94.  Giesel,  Renger,  and  Strebel  (1973)  treated  unsaturated  vertical 
flow  with  hysteresis  using  the  Crank-Nicholson  method  and  Jacobi  iteration 
to  deal  with  the  nonlinearity  of  the  equation  for  soil  moisture  movement. 

The  alternating  direction  implicit  (ADI)  method  was  used  by  Bresler  (1975) 
and  by  Busscher  (1979)  to  model  nonsteady  infiltration  from  surface  and 
subsurface  sources.  Comparison  of  the  former  model  with  laboratory  data  is 
made  by  Bresler  and  Russo  (1975).  Two-dimensional  flow  was  considered  in 
each  model,  as  implied  by  use  of  the  ADI  method.  De  Jong  and  Cameron  (1979) 
used  an  explicit  difference  formulation  of  the  diffusion  form  of  the  soil 
moisture  equation  with  explicit  linearization  of  the  diffusivity  and  con¬ 
ductivity  to  study  the  water  movement  through  soils  with  a  field  crop.  Root 
extraction  and  interception  of  precipitation  by  vegetation  were  included. 


and  soil  water  dynamics  in  layered  soils  (Hillel  and  Talpaz,  1977).  These 
four  papers  are  a  tribute  to  the  potential  of  models  for  meaningful  numerical 
experiments.  While  much  of  the  simulation  output  is  simply  reasonable, 

quantitative  comparisons  of  the  effects  of  single  parameter  variation  are 

possible.  Such  single  parameter  variation  is  essentially  impossible  in 
physical  experiments,  but  it  is  demonstrably  necessary,  as  recognized  by  most 
modelers,  to  tie  the  calculations  down  as  frequently  as  possible  to  laboratory 
and  field  data.  The  "typical  soils"  approach  of  this  series  of  papers, 
however,  may  prove  to  be  the  only  rational  way  to  treat  global  soil  moisture 
problems . 

90.  Hillel,  Van  Beek,  and  Talpaz  (1975)  discuss  the  relative  character¬ 

istics  and  merits  of  microscopic  (single  root)  and  macroscopic  (bulk)  models 
of  water  extraction  from  the  soil  by  roots  and  then  develop  a  microscopic 

model  of  the  phenomenon.  Their  CSMP  model  includes  solute  movement.  The 

model  is  formulated  in  cylindrical  coordinates,  as  were  the  analytical 
models  of  Gardner  (1960)  and  Cowan  (1965) .  Root  water  extraction  under  two 
transpiration  demand  rates  was  included.  Soil  water  potentials  were 
calculated  for  a  basic  case,  and  for  different  initial  moisture  contents, 
different  transpirational  demand,  and  soil  resistance  values.  The  effects 

of  rooting  depth  and  density  on  soil  moisture  were  simulated  with  a  macro¬ 
scopic  model  by  Hillel,  Talpaz,  and  Van  Keulen  (1976),  while  the  latter 
model  was  modified  by  Hillel  and  Talpaz  (1976)  to  include  effects  of  root 
growth  and  death.  Growth  was  treated  as  root  extension  and  root  prolifera¬ 
tion  within  a  volume.  The  model  was  compared  to  laboratory  data,  along  with 
the  model  of  Molz  and  Remson  (1970),  by  Belmans,  Feyen,  and  Hillel  (1979) 
and  by  Feyen,  Belmans,  and  Hillel  (1980).  The  results  of  the  comparison 
were  generally  satisfactory  for  both  models.  Details  are  given  in  the  cited 
papers . 

91.  Lambert  and  Penning  de  Vries  (1973)  and  Van  Bavel  and  Ahmed 
(1976)  published  models  using  CSMP  to  treat  the  entire  soil-plant-atmosphere 
system.  The  former  model  used  the  microscopic  approach,  while  the  latter 
used  the  macroscopic  approach,  to  soil  water  extraction  by  roots.  Each 
model  considers  details  of  the  canopy,  including  heat  balances  for  the 
leaves  to  better  treat  transpiration.  Comparison  of  maximum  evapotrans- 
piration  rates  for  their  simulated  sorghum  crop  with  generally  accepted 
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hydraulic  conductivity  via  constant  flux  procedures.  Calculations  were  made 
with  wetting,  drying,  and  scanning  curves  to  consider  hysteresis  in  both 
hydraulic  conductivity  and  soil  water  tension  relations  to  moisture  content. 
Correspondence  of  model  results  with  measurements  were  improved  when 
hysteresis  effects  were  considered,  as  opposed  to  single-valued  relations 
from  only  wetting  or  drying  of  the  soil. 

87.  Van  Keulen  and  Hillel  (1974)  used  a  CSMP-based  model  to  evaluate 
the  effects  of  vapor  diffusion  at  very  low  soil  moisture  contents,  applying 
the  technique  of  Hanks  and  Gardner  (1965).  Hillel  (1975)  simulated  the 
soil  moisture  and  evaporation  rate  response  to  cyclic  variation  of  potential 
evaporation  (called  "evaporativity"  by  Hillel) ,  calculating  larger  soil 
moisture  contents  relative  to  steady  potential  evaporation  after  several 
days.  Recovery  of  near-surface  soil  moisture  during  nighttime  suspension  of 
evaporation  was  modeled.  (Most  of  the  results  of  these  and  following  papers 
authored  or  coauthored  by  Hillel  are  included  in  Hillel  (1977).) 

88.  Van  Bavel  and  Hillel  (1976)  developed  a  comprehensive  treatment 
of  evaporation,  considering  both  water  and  heat  transfer  in  the  soil,  as 
well  as  energy  forcing  and  aerodynamic  transfer  in  the  near-surface  atmosp¬ 
here.  They  found  that  bare  soil  does  not  support  the  concept  of  potential 
evaporation  very  well,  as  evaporation  continued  to  fall  due  to  a  number  of 
feedback  effects  (such  as  albedo,  emissivity,  and  temperature),  the  change 
in  the  first  four  days  amounting  to  12  percent  of  the  first  day's  evapora¬ 
tion.  They  also  found  that  calculations  with  the  Van  Bavel  formula  for 
potential  evaporation  (Van  Bavel,  1966)  were  sufficiently  similar  to  model 
calculations  as  to  challenge  the  wisdom  of  the  more  complex  approach.  (If 
one  is  satisfied  with  potential  evapotranspiration,  this  is  true,  but  com¬ 
pare  Hanks  et  al.  (1973)  and  Appendix  B.)  Van  Bavel  and  Hillel  applied 
their  model  to  seven  locations  throughout  the  continental  United  States 
using  June  data,  and  discussed  variations  due  to  the  differing  diurnal 
atmospheric  conditions .  This  paper  is  recommended  reading  as  an  excellent 
discussion  of  bare  soil  evaporation. 

89.  Hillel,  Van  Bavel,  and  Talpaz  (1975)  simulated  evaporation  from  a 
soil  covered  by  a  mulch  of  hydrophobic  aggregates.  This  model  was  later 
used  without  the  surface  mulch,  but  with  three  typical  soils  labeled  sand, 
loam,  and  clay  to  evaluate  profile  water  storage  (Hillel  and  Van  Bavel, 

1976)  hysteresis  effects  with  cyclic  potential  evaporation  (Hillel,  1976), 
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84.  Van  der  Ploeg  (1974)  simulated  one-dimensional  infiltration  into 
soils  using  CSMP.  He  extolled  the  virtues  of  the  language  for  soil  scientists 
unfamiliar  with  advanced  mathematics  or  computer  programming.  He  provides  a 
clear  discussion  of  the  formulation  of  the  equations  for  computer  applica¬ 
tion,  which  is  essentially  a  calculation  of  the  vertical  moisture  flux 
divergence  with  layer-averaged  conductivities,  and  potential  gradients 
between  layers  using  Darcy's  law.  The  flux  divergence  results  in  a  change 

in  moisture  content  in  each  layer.  He  finds  good  correspondence  between 
model  output  and  the  calculations  of  Philip  (1957a)  and  Parlange  (1971)  for 
Yolo  light  clay.  Van  der  Ploeg  and  Benecke  (1974)  simulated  one-,  two-, 
and  three-dimensional  infiltration  in  soils.  They  used  an  alternating- 
direction  procedure  (explicit)  for  the  two-dimensional  case,  while  the  three- 
dimensional  case  was  reduced  to  one  dimension  (radial) .  They  found  good 
correspondence  with  results  by  Philip  and  Parlange,  and  noted  that  their 
model  represented  field  data  as  well  as  the  model  of  Bresler  et  al.  (1971). 

85.  Beese,  Van  der  Ploeg,  and  Richter  (1977)  tested  their  model 
against  field  data,  a  218-day  experiment  on  fallow  loess  soil.  Field  data 
included  tensiometer  measurements  at  11  depths  to  2  m  (6  ft)  in  the  field 
and  in  a  soil  monolith,  which  had  been  removed  from  the  field  and  trans¬ 
ported  to  a  lysimeter  station  a  few  kilometers  (miles)  away.  Evaporation 
was  estimated  from  a  correlation  of  potential  evaporation/actual  evaporation 
(lysimeter)  and  matric  potential  at  5-cm  (2-in.)  depth.  Precipitation  was 
measured  daily.  The  soil  capillary  conductivity  was  measured  as  a  function 
of  depth.  Field  data  was  plotted  on  graphs  from  a  CSMP-based  model  output 
for  depths  of  20,  40,  60,  100,  and  140  cm  (8,  16,  24,  40,  and  56  in.). 

Model  accuracy  increases  with  depth,  while  the  average  departures  of 
calculated  from  measured  values  for  all  depths  were  less  than  15  percent 
(matric  potential  units) .  The  departure  may  have  been  due  to  the  values  of 
evapotranspiration  and  capillary  conductivity  determined  as  inputs,  as  well 
as  model  errors,  according  to  the  authors. 

86.  Dane  and  Wierenga  (1975)  modeled  the  effect  of  hysteresis  on 
infiltration,  redistribution,  and  drainage  in  a  layered  soil  with  CSMP. 

They  compared  model  output  with  laboratory  data  for  moisture  contents  in 
Glendale  clay  loam  over  river  sand.  Soil  water  tension  was  measured  via 
tensiometers,  soil  moisture  content  via  neutron  scatter  techniques,  and 
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This  equation  was  also  solved  via  the  Douglas-Jones  method.  The  equations 
were  linked  via  use  of  the  water  table  location  as  a  lower  boundary  condition 
on  the  unsaturated  zone,  and  defintion  of  S  as  soil  water  in  the  unsaturated 
zone  which  was  unavailable  for  evaporative  loss.  Iteration  was  not  necessary 
for  small  time  steps  when  the  saturated  zone  equation  was  solved  for  the 
water  table  height  using  values  of  S  and  w  from  the  prior  time  step,  but 
it  could  be  used  when  necessary  or  to  achieve  larger  time  increments. 

80.  Comments  on  the  Pikul,  Street,  and  Remson  paper  were  published  by 
Vachaud  and  Vauclin  (1975)  with  a  reply  by  Pikul,  Street,  and  Remson  (1975). 

In  comparison  to  laboratory  experiment  data,  the  former  authors  found  larger 
horizontal  water  movement  in  the  unsaturated  zone  than  allowed  by  the  assump¬ 
tions  of  the  model  and  strongly  challenged  the  concept  of  the  storage  coeffi¬ 
cient,  S  .  The  reply  was  essentially  “application  to  a  watershed  cannot  be 
made  with  the  same  approach  used  in  the  laboratory  because  of  sheer  size." 

S  was  claimed  to  be  useful,  while  lateral  flow  in  the  unsaturated  zone 
would  likely  be  far  less  than  in  the  laboratory  experiment. 

81.  A  number  of  models  have  been  formulated  using  the  IBM  CSMP,  which 
was  developed  for  use  by  individuals  generally  unfamiliar  with  computer 
programming  (IBM  Corporation,  1972) .  These  models  are  reviewed  as  a  group 
in  the  following  paragraphs,  as  the  papers  share  a  common  modeling  techni¬ 
que,  albeit  developed  by  IBM  rather  than  originated  by  an  individual  soil 
scientist  with  programming  skills. 

82.  All  CSMP  models  discussed  in  the  following  paragraphs  use 
explicit  integration  procedures,  using  either  the  Milne  method  or  a  fourth- 
order  Runge-Kutta  scheme.  The  user  of  these  methods  is  protected  from 
instability  in  the  explicit  integrations  by  software  checks,  including 
reference  to  a  user-specified  accuracy  requirement.  Complete  programs  have 
been  published  for  most  of  the  models,  as  they  are  quite  short  because  of 
the  amount  of  detail  left  up  to  the  software. 

83.  Van  Keulen  and  Van  Beek  (1971)  simulated  water  movement  in 
layered  soils  with  CSMP  to  evaluate  the  influence  of  a  plow  layer  and 
tillage  hardpan.  They  included  appendices  which  discuss  the  influence  of 
layer  thickness  selected  for  the  vertical  model  resolution  and  the  magnitude 
of  the  time  step  used  in  integration.  Bhuiyan  et  ai.  (1971)  used  CSMP  to 
model  vertical  infiltration  into  unsaturated  soils.  They  found  fairly  good 
agreement  between  their  model  and  the  power  series  solution  due  to  Philip 
(1957a)  for  Yolo  light  clay. 
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This  problem  is  solved  with  both  a  Taylor  series  expansion  and  the  method  o£ 
Douglas  and  Jones,  described  above.  They  found  that  the  Taylor  series 
required  too  high  an  accuracy  in  higher  order  derivatives  of  the  soil 
moisture  characteristics  curves,  which  are  difficult  to  determine,  and  that 
it  was  not  acceptable  for  this  reason.  Solution  with  the  Douglas- Jones 
predictor-corrector  method  of  finite  differencing  was  satisfactory  in 
comparison  with  experimental  data  from  gravity  drainage  of  a  column  of  soil. 

77.  Pikul,  Street,  and  Remson  (1974)  published  a  numerical  model  for 
two-dimensional  applications  based  on  coupled  one-dimensional  models  for  the 
unsaturated  and  saturated  zones.  While  known  limitations  to  this  approach 
were  cited,  the  potential  benefits  of  a  relatively  efficient  application  to 
watershed-size  problems  were  cited  as  justification.  Further,  at  the  water¬ 
shed  scale,  many  of  the  assumptions  are  reliable. 

78.  The  equation  of  soil  moisture  flow  in  the  unsaturated  zone  was 
written  in  terms  of  pressure  head  using  the  specific  moisture  capacity, 

C  ,  as 


where 

z  =  vertical  coordinate 

K  =  hydraulic  conductivity 

¥  *  pressure  head  (matric  potential) 

C  =  specific  moisture  capacity  (=80/3'}') 
t  =  time 

This  equation  was  solved  using  the  Douglas-Jones  method  discussed  in  the 
earlier  papers. 

79.  The  saturated  flow  zone  was  modeled  with  the  equation 

Ko  -  fe  (h  ■  s<*-«  It  -  (CU) 

where 

Kq  =  saturated  hydraulic  conductivity 
x  =  horizontal  independent  variable 
h  =  height  of  the  water  table  above  datum 

S  =  a  storage  coefficient  used  to  link  the  saturated  and  unsaturated 
zones 

t  =  temporal  independent  variable 
w  =  a  sink  or  source  term 
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73.  The  extraction  term  was  first  written  as  a  function  of  only  depth 
and  transpiration  rate.  The  resulting  equation  was  solved  with  an  Adams 
predictor-corrector  method  started  by  a  Runge-Kutta  procedure.  Reasonable 
agreement  was  found  with  measured  moisture  flux.  A  second,  more  realistic 
model  with  the  extraction  term  as  a  function  of  moisture  content,  as  well  as 
transpiration  rate  and  depth  was  then  formulated.  Expressions  for  both  the 
rooting  density  function  and  the  resulting  root  extraction  function  with 
depth  are  given  in  the  cited  reference.  The  Douglas- Jones  predictor- 
corrector  method  was  used  to  solve  the  resulting  mathematical  model.  Good 
results  were  found  for  the  first  several  days  of  simulation  in  comparison  to 
data  from  a  field  experiment,  but  the  simulation  and  field  data  departed 
strongly  after  that  period.  There  was  evidence  in  the  data  that  maximum 
extraction  was  shifting  to  depths  of  lower  root  uptake  capacity  as  the  soil 
dried,  a  feature  not  included  in  the  model  (see  Belmans,  Feyen,  and  Hillel 
(1979)  for  additional  comparison  to  data) . 

74.  The  Douglas- Jones  predictor-corrector  method  is  outlined  in  an 
Appendix  to  the  paper  (see  also  Remson,  Hornberger,  and  Molz,  1971).  In  the 
Douglas-Jones  method  a  predictor  equation  is  used  to  evaluate  the  unknowns 
at  one-half  step  forward  in  time  using  an  implicit  formulation  and  solution 
via  the  tridiagonal  matrix  method.  These  values  are  then  used  to  derive  the 
coefficients  at  the  half-time  step  for  use  in  a  corrector  equation.  The 
latter  is  again  implicit,  but  it  is  applied  over  the  entire  time  step. 
Solution  is  again  by  the  tridiagonal  matrix  method.  The  nonlinearity  of  the 
equations  during  the  corrector  step  is  therefore  removed  by  use  of  known 
values  from  the  predictor  equation  solution,  while  for  the  predictor  step 
the  coefficients  are  taken  as  appropriate  to  the  beginning  of  the  step 
(explicitly).  First  derivatives  are  written  explicitly  for  each  step,  also. 
The  method  is  convenient,  efficient,  and  accurate. 

75.  Molz  and  Remson  (1971)  apply  the  1970  model  to  a  number  of  cases 
of  soil  moisture  extraction  by  roots,  but  unfortunately  without  adequate 
comparison  to  experimental  data. 

76.  Hornberger  and  Remson  (1970)  model  the  one-dimensional  saturated- 
unsaturated  transient  flow  problem  with  a  formulation  of  the  basic  equation 
in  pressure  head.  They  depart  from  the  approaches  of  Rubin  (1968)  and 
Freeze  (1969)  discussed  above,  however,  in  that  they  model  a  discontinuity 
in  pressure  head  at  the  water  table.  Their  rationale  is  given  in  the  paper. 
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70.  Remson,  Fungaroli,  and  Hornberger  (1967)  present  still  another 
approach  to  numerical  modeling  of  soil  moisture.  They  formulated  a  three- 
dimensional  model  based  on  flow  divergence,  V  V  =  30/3t  ,  in  which  the 
six  components  of  low  toward  and  away  from  a  point  in  three  dimensions  were 
first  written  out  in  terms  of  head,  conductivity,  and  flow  area.  Because 
discrete  differences  between  adjacent  spatial  grid  points  were  used  in 
writing  the  volumetric  flow  rate  components,  the  finite  difference  approxima 
tion  was  complete  with  the  addition  of  an  explicit,  backward-difference 
term  for  the  time  rate  of  change  of  volumetric  moisture  content.  Tabulated 
values  of  K  versus  6  and  h  versus  0  were  used  with  an  iteration 
scheme  to  compute  the  solution,  i.e.,  the  spatial  variation  of  0  at  each 
time  step. 

71.  The  model  was  applied  to  a  problem  of  a  draining  soil  with 
evaporation,  but  in  only  one  dimension.  Comparison  with  computed  results 
for  a  similar  case  by  Remson  et  al.  (1965)  was  good  without  hysteresis,  but 
computed  results  with  hysteresis  had  a  problem  near  the  region  of  transition 
between  wetting  and  drying,  which  was  attributed  to  use  of  an  insufficient 
number  of  hysteretic  scanning  curves  between  the  wetting  and  drying  curves. 

72.  Molz  and  Remson  (1970)  formulated  a  mathematical  model  for  soil 
water  extraction  by  roots.  Contrary  to  the  approach  of  an  earlier  paper — 
Molz  et  al.  (1968) — in  which  flow  to  a  single  root  was  modeled  (microscopic 
model),  they  used  the  macroscopic  approach.  In  this  case,  the  bulk  char¬ 
acteristics  of  the  roots  and  their  extraction  of  water  with  depth  are 
treated.  A  negative  source  term  was  added  to  the  equation  of  continuity, 
which  was  then  combined  with  Darcy's  law  to  derive  an  equivalent  of  the 
Richards  equation,  but  with  the  added  source.  The  diffusivity  modification 
due  to  Klute  was  then  used  to  derive  the  final  equation: 

!£  =  -1  (D  ii)  _  M  _  s  (cl2 
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where 


0  =  volumetric  moisture  content 
t  =  time 

z  =  the  vertical  coordinate  (positive  down) 
D  =  soil  moisture  diffusivity 
K  =  soil  moisture  conductivity 
S  *  the  root  extraction  source  term 


After  an  informative  discussion  of  the  problem  of  incorporating  treatment  of 
both  saturated  and  unsaturated  domains  in  a  single  model,  and  of  the  physics 
of  ground  water  recharge  and  discharge  relative  to  water  table  depths. 

Freeze  formulates  his  model  in  terms  of  pressure  head  (matric  potential) 
using  a  Crank-Nicholcon  implicit  method  with  linearization  of  coefficients 
via  extrapolation  after  Rubin  and  Steinhardt  (1963) .  The  system  of  simult¬ 
aneous  equations  is  solved  by  the  tridi. g'nal  matrix  approach  frequently 
used  (Richtmyer,  1957). 

67.  Freeze  (1969)  applied  the  model  to  several  hypothetical  situations 
using  experimental  data  for  characteristic  curves  of  three  soils.  A  sub¬ 
sequent  paper  was  noted  in  which  comparison  of  results  with  field  data  would 
be  made,  but  it  was  not  available  for  review.  Nonetheless,  this  paper  is 
highly  recommended  reading  for  its  development  of  the  problem. 

68.  Schneider  and  Luthin  (1978)  utilized  the  Rubin  (1968)  formula¬ 
tion  of  the  saturated-unsaturated  flow  problem  to  study  perched  water 
tables.  Their  problem  required  addition  of  a  source  term,  which  was 
included  in  the  finite  difference  equations  implicitly.  The  solution  was 
obtained  by  iteration  as  used  by  Rubin.  Model  results  were  compared 
with  data  from  a  laboratory  experiment.  The  results  indicated  that  the 
model  was  calculating  too  large  an  unsaturated  horizontal  flow,  but 
qualitative  comparison  of  both  perched  and  ground-water  tables  was 
valuable. 

69.  In  an  earlier  paper,  Taylor  and  Luthin  (1969)  had  approached 
the  combined  unsaturated-saturated  zone  problem  by  first  solving  for  the 
two  zones  separately,  then  adjusting  the  boundary  condition  beteen  the 
zones.  Because  they  applied  their  numerical  model  to  drawdown  of  an 
aquifer  by  a  well,  the  equations  were  formulated  in  cylindrical  coordin¬ 
ates.  Expressions  for  both  the  saturated  and  unsaturated  zones  were 
written  in  terms  of  hydraulic  head  $  =  h  +  z  (=H)  ,  with  the  difference 
being  the  absence  of  a  time  derivative  in  the  saturated  zone.  Explicit 
finite  difference  formulation  was  used  and  solved  by  iteration.  Variable 
spacing  was  used  in  the  vertical  and  radial  directions  to  balance  computa¬ 
tional  accuracy  and  efficiency.  Good  correspondence  was  found  with 
laboratory  data  for  sand. 
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62.  The  alternating-direction  implicit  (ADI)  method  of  Peaceman  and 
Rachford  (1955)  was  used  for  the  solution  of  the  two-dimensional  infiltra¬ 
tion  problem.  In  this  method  (see  also  Remson,  Hornberger,  and  Molz,  1971) 
the  two-dimensional  equation  is  solved  in  two  steps  per  advance  in  time, 
which  may  be  either  one-half  time  increments  or  whole  time  increments.  In 
the  latter  case,  two  time  increments  are  required  for  complete  application 
of  the  method.  Basically,  the  approach  is  to  first  write  the  equation 
implicitly  in  one  variable,  say  x  ,  with  the  variation  of  the  second  vari¬ 
able,  say  z  ,  denoted  explicitly.  This  results  in  a  system  of  equations 
for  H  at  x-grid  points  which  is  solved  via  tridiagonal  matrix  techniques. 
The  second  time  step  is  then  accomplished  by  reversing  the  implicit-explicit 
application  to  the  variables,  so  that  x  variations  are  written  explicitly 
and  z  variations  implicitly.  The  new  system  of  equations  is  again  solved 
via  a  tridiagonal  technique  for  the  z-grid  points,  and  one  cycle  is  complete 
This  is  an  efficient,  stable  technique  for  solution  of  problems  in  two  space 
dimensions . 

63.  Rubin  (1968)  uses  the  extrapolation  technique  of  Rubin  and 
Steinhardt  (1963)  to  linearize  the  equations  by  explicit  calculation  of  the 
Y (v)  and  Z(v)  coefficients  obtained  with  the  Kirchhoff  transformation. 

64.  The  falling-water-table,  ditch-drainage  problem  is  also  treated 

by  an  ADI  scheme,  but  because  the  saturated  region  at  the  base  results  in  an 

elliptic  equation  (no  time  variation),  iteration  is  required  for  a  solution. 

An  iteration  term  of  the  form  I  KAH  is  added  to  the  left  side  of  the 

m 

finite  difference  forms  (x-implicit  and  z-implicit)  of  Equation  Cll  and  the 

calculations  are  repeated  for  each  complete  time  step  until  two  consecutive 

iterations  are  sufficiently  close.  The  iteration  parameter  is  cycled  during 

iteration  to  improve  convergence.  It  is  given  by  I  =  Rs  ,  where  s  =  0, 

m 

1,  2,  ...  S  ,  and  R  and  S  are  constants.  The  values  used  in  Rubin 
(1968)  were  R  =  0.22  and  S  =  6  .  The  cited  reference  should  be  consulted 
for  details  and  complete  expressions  of  the  finite  difference  equations. 

65.  The  results  of  this  model  were  again  found  to  be  reasonable  and 
enlightening,  but  experimental  data  was  not  available  for  evaluation. 

66.  Freeze  (1969)  used  a  model  similar  to  that  of  Rubin  and 
Steinhardt  (1963)  to  study  one-dimensional,  vertical,  unsteady,  unsaturated 
flow  above  a  recharging  or  discharging  ground-water  flow  system.  His  table 
of  prior  work  and  discussion  thereof  is  frequently  cited  by  later  authors. 


the  method  of  computation  cannot  change  the  process  physics.  Details  are 
marked,  and  frequently  merged,  by  use  of  correlation  relations  for  the  unknown 
physics  of  the  processes  of  moisture  increase  and  decrease.  The  process  is 
nonetheless  complex,  and  simple  models  are  simply  wrong  in  cases  which 
depart  in  too  large  measures  from  the  data  used  to  define  the  model  relations. 

106.  On  the  other  hand,  despite  great  strides  during  the  past  several 
decades,  many  important  elements  of  the  physics  of  moisture  exchange 
processes  remain  unknown,  or  require  inaccessible  detail  or  extensive 
computer  time  to  treat.  Thus,  the  models  of  the  foregoing  section  may  not 
produce  superior  results  in  comparison  to  a  simple  budget  model  in  real 
cases  which  are  bounded  by  temporal  and  financial  resources.  Except  for  the 
quantity  of  data  required  in  some  applications,  the  budget  models  may  be 
implemented  with  far  less  computing  power  than  the  physically  more  accurate 
numerical  models. 

107.  Thornthwaite  and  Mather  (1954)  published  an  updated  version  of  a 
budget  (mass  balance)  model  for  soil  moisture  prediction  with  application  to 
soil  tract ionability .  Water  is  input  to  the  soil  column  via  measured 
precipitation,  while  it  is  removed  by  both  gravitational  drainage  and 
evapo transpiration  at  moisture  contents  above  field  capacity,  and  by  only 
evapotranspiration  below  this  amount.  The  Thornthwaite  method  of  calcula¬ 
ting  potential  evapotranspiration  is  used  (see  Appendix  B) ,  while  actual 
evapotranspiration  is  adjusted  to  reflect  the  soil  moisture  remaining  via  a 
linear  relationship.  Soil  type  is  considered  in  determination  of  field 
capacity,  but  other  properties  such  as  conductivity  are  not  used.  They 
claim  accurate  determination  of  moisture  content  and  flow  through  soils,  and 
compare  to  field  data. 

108.  A  budget  model  developed  by  personnel  of  the  US  Forest  Service 
and  US  Army  Corps  of  Engineer  Waterways  Experiment  Station  is  discussed  by 
Burke  and  Turnbull  (1959) .  This  model  is  discussed  more  fully  in  the  body 
of  this  report,  but,  briefly,  it  utilizes  precipitation  as  a  model  input 
which  is  distributed  into  two  15-cm  (6-in.)  layers  as  determined  by 
statistically  derived  accretion  relations,  while  depletion  of  the  layers 
is  effected  by  moisture-content-dependent  depletion  relations,  also 
statistically  derived  from  a  large  data  set.  The  relations  are  different 

for  sand,  silt,  and  clay  soils,  and  for  winter,  transition,  and  summer  seasons. 


The  model  avoids  the  questionable  concept  of  field  capacity  by  use  of  a  field 
maximum  water  content,  which  is  again  derived  from  the  data.  Good  results 
were  obtained  for  well  drained,  fine-grained  soils  in  locales  similar  to  those 
of  the  original  data;  however,  poorer  correspondence  was  obtained  between 
measurement  and  model  calculations  in  poorly  drained  soils,  soils  with  high 
organic  content,  and  tropical  regions. 

109.  Baier  and  Robertson  (1965  and  1966)  developed  a  six- layer  soil 
moisture  budget  model  in  which  primary  depletion  was  effected  via  evapo- 
transpiration.  Evapotranspiration  was  determined  by  means  of  correlations 
using  up  to  six  input  variables,  namely,  maximum  temperature,  temperature 
range,  solar  energy  at  the  top  of  the  atmosphere,  total  solar  energy  at  the 
surface,  wind  run,  and  vapor  deficit.  Depletion  from  each  level  was  allowed 
through  a  factor  accounting  for  soil  and  plant  root  characteristics,  and  a 
linear  relation  between  actual  evapotranspiration  and  layer  moisture  content, 
as  in  the  Thornthwaite  approach.  Runoff  was  estimated  via  a  relation  which 
included  rainfall  and  soil  moisture  in  the  upper  layer.  Daily  evapotrans¬ 
piration  was  calculated  for  days  with  rain  before  adding  the  appropriate 
precipitation  amount.  Baier  (1971)  found  the  Baier  and  Robertson  Versatile 
Soil  Moisture  Budget  (VSMB)  to  produce  results  closer  to  the  output  of  the 
Penman  equation  for  potential  evapotranspiration  than  did  the  Thornthwaite 
approach.  The  VSMB  has  been  used  in  several  other  models.  Hildreth  (1977) 
provides  a  discussion  and  error  analysis  of  this  model  in  connection  with 
outlines  of  several  other  mass  balance  approaches.  Hildreth  also  outlines 
several  models  for  evapotranspiration  and  potential  evapotranspiration,  also 
discussed  in  Appendix  B. 

110.  A  relatively  simple  budget  model  for  irrig  tion  scheduling  was 
published  by  Jensen,  Robb,  and  Franzoy  (1970).  Essentially  only  rainfall 
and  evapotranspiration  are  considered,  although  the  Penman  approach  with 
active  consideration  of  varying  crop  coefficients  through  the  growing  season 
was  used.  Both  prior  weather  and  forecasts  were  used  over  a  period  of 
several  days  in  model  applications  to  estimate  moisture  deficits  and 
required  irrigation. 

111.  A  combination  budget-physical  process  model  was  developed  by 
Jones  and  Verma  (1971).  Each  model  day,  rainfall  minus  potential  evapora¬ 
tion  was  simply  distributed  into  a  layered  soil,  bringing  each  layer  to 
saturation  in  turn,  until  the  water  was  exhausted.  The  soil  moisture  was 


then  linearly  redistributed  into  the  soil  layers  while  holding  either  the 
lowest  layer  moisture  content  constant  or  the  surface  layer  at  saturation, 
depending  on  the  depth  of  penetration  of  the  wetting  front.  As  long  as  the 
surface  layer  remained  above  air-dry  moisture  content,  all  evaporation 
moisture  was  taken  from  the  uppermost  layer  and  was  set  equal  to  potential 
evaporation  (determined  as  0.7  times  evaporation  measured  by  US  Weather 
Bureau  Class  A  evaporation  pan) .  When  the  surface  moisture  content  fell 
below  the  air-dry  value,  evaporation  was  calculated  with  the  general  approach 
of  Hanks  and  Bowers  (1962) ,  but  applied  to  the  flow  equation  written  in 
terms  of  moisture  content  and  moisture  diffusivity.  Good  correspondence  was 
found  between  model  simulation  and  field  data.  The  importance  of  accuracy  in 
diffusivity  values  near  saturation  was  noted,  in  keeping  with  prior  work 
noted  above.  Predicted  values  of  soil  moisture  content  were  generally 
within  10  percent  of  measurements  over  a  period  of  A3  days  of  natural 
weather . 

112.  Stuff  and  Dale  (1978)  used  a  budget  model  to  account  for  shallow 
water  table  influences  on  soil  moisture  under  corn.  Actual  evapotranspira- 
tion  was  estimated  from  Class  A  pan  data  which  was  adjusted  for  both  a  crop 
factor  dependent  on  stage  of  growth  and  nonmoisture  stress  factors  related 
to  demand  rate  and  moisture  deficits.  Empirical  relations  were  developed 
for  water  table  depth  and  capillary  rise  from  two  years  of  field  data  and 
then  utilized  to  test  the  resulting  model  against  two  different  years'  data. 
The  general  performance  of  the  model  was  considered  adequate,  with  concern 
expressed  about  the  validity  of  the  data  used  in  deriving  the  empirical 
relationships  used,  including  the  range  of  moisture  contents  represented  in 
the  original  data  set. 

113.  Lettau  (1969)  divided  precipitation  input  to  a  surface  into  two 
main  components  in  the  formulation  of  a  budget-process  hybrid  model.  That 
portion  of  the  precipitation  which  ran  off  or  was  evaporated  from  the 
surface  during  the  input  data  interval  was  parameterized  in  a  budget  model, 
and  subtracted  from  the  total  precipitation,  resulting  in  a  reduced  forcing 
function  for  trends  of  soil  moisture.  The  balance  of  actual  runoff  and 
evapotranspiration  for  a  given  data  interval  was  assumed  to  come  from  water 
which  had  been  stored  in  the  soil  during  previous  intervals.  The  immediate 
fluxes  were  modeled  in  terms  of  the  dual  forcing  functions  of  precipitation 


and  solar  energy,  while  the  delayed  fluxes  were  set  proportional  to  soil 
moisture  content.  This  approach  resulted  in  a  simple  differential  equation 
for  the  trends  of  soil  moisture,  which  was  solved  by  an  unspecified  numerical 
integration  technique.  The  model  and  parameterization  were  intended  for 
watershed  application  over  periods  of  several  months,  and  results  for  the 
North  American  Midwest  were  found  to  be  superior  to  other  large-scale  models. 
The  model  approach  was  used  by  Lettau  and  Baradas  (1973)  for  a  small  water¬ 
shed  in  the  Phillipine  Islands,  by  K.  Lettau  (1974)  in  a  steppe  environment, 
by  Lettau  and  Lettau  (1975)  for  the  tundra  and  boreal  forests  of  Canada,  and 
by  Hall  (1977)  for  the  North  American  Great  Plains.  In  the  latter  two 
applications,  the  evapotranspiration  model  was  used  in  conjunction  with 
solar  irradiance  and  surface  energy  budget  models  to  address  both  input  and 
distribution  of  energy  at  the  surface  of  the  earth. 

114.  Each  of  the  mass  balance  approaches  has  had  to  rely  on  para¬ 
meterization  or  empirical  relations  to  incorporate  process  physics  into  the 
model  equations.  An  attempt  to  avoid  this  was  made  by  Lund  and  Needleman 
(1974) ,  who  developed  straightforward  regression  relations  between  meteor¬ 
ological  data  and  measured  soil  moisture  in  the  layer  from  15  to  30  cm  (6  to 
12  in.)  for  a  tropical  location.  Even  with  consideration  of  time  lags,  they 
were  only  able  to  distinguish  wet  and  dry  seasons,  with  the  possibility  of 
determining  rates  of  transition  noted.  It  is  of  value  to  eliminate  an 
unworkable  method,  as  these  authors  have  done,  so  effort  can  be  directed  to 
less  simple  methods  without  hope  for  a  raw  statistical  approach  continually 
surfacing.  Soil  moisture  modeling  has  proven  to  be  a  strongly  nonlinear 
problem  of  much  greater  complexity  than  one  which  could  be  solved  so  simply. 
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